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The Solution of Aeroelastic Problems by 
Electronic Analogue Computation’ 


JONATHAN WINSON?t 


Fairchild Engine and Airplane Corporation 


SUMMARY 


The ease with which complex physical problems are resolved 
by machine computation of the analogue type has led to an in- 
spection of flutter and related aeroelastic phenomena with a view 
toward automatic solution of these problems. It is found that 
automatic solution is feasible and can be carried out in the follow- 
ing manner: The basic aeroelastic equations for arbitrary wing 
motion are written and entered in a commercial electronic com- 
puter. The elastic, inertia, and aerodynamic properties of the 
physical system are translated into fixed potentiometer settings 
of the computer, a small number of variable settings being re 
served for aircraft forward velocity and initial flight conditions. 
The effects of the vortices in the wake of the wing are represented 
by the outputs of simple electrical filters. For the understand 
ing of the flutter characteristics of an aircraft, say, machine runs 
are made at several forward velocities, the wing in each run being 
disturbed initially. Pen traces of the resulting wing oscillations 
as functions of time yield all required dynamic information. The 
oscillations diverge above the flutter speed. To illustrate the 
application, sample results with appropriate numerical com 
parisons are presented. 

The analogue computation technique recommends itself by 
its flexibility and by the large timesaving achieved through its 
use. Systems with many degrees of freedom and systems con- 
taining nonlinear components, both difficult to treat mathe- 
matically, can be analyzed by machine with comparative ease 


INTRODUCTION 


Ov MAY SAY THAT the purpose of the study of air- 
craft dynamics is to supply to the aircraft de- 
signer, at an early stage of the design process, detailed 
information as to the flight characteristics and strength 
The fact that this 
information is needed as the design is conceived, long 
before construction and flight, requires that the study, 


requirements of his proposed craft. 


Presented at the Aeroelasticity Session, Eighteenth Annual 
Meeting, I.A.S., New York, January 23-26, 1950. 

*The author wishes to express his thanks to S. Fifer, H. H. 
Adise, and others of the Reeves Instrument Company and to 
H. Levenstein, of Fairchild, for their valuable aid in carrying 
out this work . 

} In Charge of Dynamics, Guided Missiles Division. 


in general, be analytic. Though close analysis is not 
synonymous with tedious effort, the aircraft dynami- 
cist might well consider it so in the light of his practical 
experience. The nature of the aerodynamics of the air- 
plane is, indeed, such as to make analysis a compli- 
cated matter. 
frequency oscillations of aerodynamic surfaces are in 


This is particularly true where high- 


volved and where consideration must be given to the 
effects of the vortices in the wakes of the surfaces. 

It has been the practice to treat the general dynamic 
problem of the action of the elastic air frame in flight 
as two problems: one, of rigid body motion; the second, 
of elastic motion of certain lifting surfaces with re- 
stricted rigid body interaction. The problem of rigid 
body motion is termed the stability and control prob- 
lem. Its treatment is usually simplified by the assump- 
tion that the previously mentioned wake effects are 
unimportant. The second problem is the aeroelastic 
one. The wake effects are considered here and com- 
plicate the mathematical treatment to a point where the 
analyst is frequently content merely to explore the sys- 
tem under consideration for a point of neutral stability 
(a flutter point), although a deeper understanding of 
the aeroelastic properties of the system is desired. 

Recent trends in design complicate this situation still 
further, making more and more untenable the assump- 
tions that allow separation of rigid body and elastic 
body effects. The field 
to the other, and a more or less unified approach to the 


work of each draws closer 
total dynamic problem, with its formidable mathe- 
matical complexity, may soon be required. 

It is the purpose of this paper to present for considera 
tion a new method of solution of the dynamic problems 
of aircraft—one that shows promise, in preliminary in- 
vestigation, of reducing considerably the difficulties of 
analysis, with little if any sacrifice in accuracy. The 
method is that of analysis by electronic analogue com- 


putation. As might be inferred from the name, it con- 


385 








386 JOURNAL OF THE 
sists of representing the physical aerodynamic and 
elastic system by another system, made up of anal- 
ogous electrical components, and of determining the 
action of the original system by observing the char- 
acteristics of the analogous one. The translation from 
one system to the other is carried out in an orderly, 
rigorous manner, common mathematical description 
of the two serving as a link between them. A com- 
mercial analogue computer, the REAC (Reeves Elec- 
tronic Analogue Computer), is used to set up the anal- 
ogous system. The special aeroelastic problem of 
flutter is treated here to illustrate the technique and 
evaluate the worth of the method. Flutter, with its 
complete consideration of the wake effects, contains the 
kernel of all of the more general dynamic problems. 

In Section I below, the aeroelastic equations of mo- 
tion for a two degree of freedom system are written in 
a form suitable for analogue computation. Two- and 
three-dimensional cases are considered, both cases with- 
out structural damping and the three-dimensional case 
with an untapered wing. Section II gives some com- 
ment on the development of analogue computation, a 
description of the REAC, and the analogue schemes cor- 
responding to the previously derived equations. Sec- 
tion III presents the machine results with numerical 
comparisons. Section IV discusses the application of 
the method to broader problems and comments on ana- 
logue setups for systems including tapered wings, struc- 
tural damping, and additional degrees of freedom. 


MATHEMATICAL STATEMENT OF THE 
PROBLEM 


SECTION I 


NOMENCLATURE 


Dimensional Quantities 


U = forward velocity of the aircraft, ft. per sec. 
p = density of the air, slugs per cu.ft. 
b, = reference semichord of the wing, ft. 


Dimensionless Quantities 


b = semichord of the wing at any section 

a = distance of elastic center of any section aft of 
mid-chord as a fraction of semichord 

x = spanwise coordinate of section measured from 


centerline of aircraft 


l = semispan of wing 

h = displacement of elastic center from reference 
plane, positive downward 

0 = angular displacement of section from refer- 
ence direction, positive nose up 

h’ = displacement of three-quarter chord point from 
reference plane, positive downward 

c = distance of center of gravity of any section aft 
of mid-chord as a fraction of semichord 

r = radius of gyration of a section about the elastic 
center 

p’ = mass of a section 

hr = parameter describing spanwise variation of h 

yo = parameter describing spanwise variation of 6 

T = time 

L = total lift on a section, positive up 

M = total moment ona section, positive nose up 


L,, Le, L; = components of lift 
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M, = component of moment 

a = kinetic energy 

V = potential energy 

Qn, Qa = generalized forces 

w = frequency of oscillation 

Wh = uncoupled bending natural frequency 

we = uncoupled torsional natural frequency 

M1, Az, Az = characteristic spanwise integrals 

m = integrated mass characteristic 

1 = integrated value of moment of inertia aboy 
the elastic center 

S = integrated value of mass static moment aboy 
the elastic center 

a}, dz, Bi, By = parameters describing lift deficiency 


Nondimensional Scheme 


The operational approach to the investigation oj 
flutter which is used may be summarized as follows 
The aeroelastic system is represented in analogue form, 
the terms in the equations of motion appearing on the 
computer as potentiometer settings. As is customary in 
three-dimensional analytic flutter work, the uncoupled 
modes of wing vibration are assumed to be known from 
other sources. At each of several forward velocities 
initial disturbances are imposed, and the reaction of the 
system is noted. The major variable from one machine 
run to the next is forward velocity, and it is clear that, 
for operation to be satisfactory from a practical point 
of view, the effect of velocity must be confined to an 
extremely small number of potentiometers. The ve 
locity, then, must not appear explicitly in more than a 
few terms in the equations of motion. This can be 
accomplished by a transformation of time coordinate 
such that variables are referred to the basic physical 
parameter, translation of the wing, rather than to real 
time. A number of schemes can be used, one of the 
most common and convenient being the adoption of 4 
nondimensional time in which real time is referred to the 
time required for translation of the wing through a dis 
tance equal toitssemichord. A transformation method 
suggested by Goland and Luke! embraces this point 
and is used here. Portions of the nomenclature and 
derivation of this section also parallel the treatment i 
reference 1. 

In the nondimensional scheme adopted, a dimen 
sionless time 7 is set equal to dimensional time ¢ divided 
by the time required for wing translation through 4 
reference semichord length. Thus, 7 = (l’/b,, where 
U is the forward velocity of the aircraft and ), 184 
reference semichord. Physical lengths are made di 
mensionless by division by 6,; forces, by division by 
pU*b,? (where p is the density of the atmosphere); fre 
quencies, by division by U/b,, ete. The result of this 
operation is to confine velocity effects to two terms tf 
the equations of motion—namely, w, and w,, the not 
dimensional natural frequencies of the fundamental 
uncoupled bending and torsional modes of the wing 
Thus, changes in two potentiometer settings of the ama: 


logue computer allow the shift from one velocity 
another. The details of this procedure are discussed 1 
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SOLUTION OF PROBLEMS 


Nondimensional notation is used in the 


Section II. 
work that follows. 


Aerodynamics 

Fig. | portrays a wing section in displaced position. 
The lift components L, Lz, and L; and the moment com- 
ponent M, are shown acting on the section. They rep- 
resent the entire aerodynamic force system. The dis 
placements of the elastic center and the three-quarter 
chord point of the section from a reference plane are 
denoted by /# and h’, respectively. The angular dis- 
placement of the section from a reference direction is 
denoted by #6. The quantity 4 is the semichord; ad and 
ch are the distances of the elastic center and center of 
gravity of the section aft of the mid-chord. 

L,, the circulatory component of the lift, acts at the 
L, and L;, apparent mass forces, 


‘ 


quarter chord point. 


dh’ 
L, = 2nb (0 = 
d 


T 


7/0 


Ly = b?[(d0/dr) + (d*h'/dr*)| (2) 
L3 = (a/2)b* (d°6/dr?) (3) 
M, = (2/2)b*(d0/dr) (4) 


In these equations, @, ds, 6;, and PB» are “‘lift defi- 
ciency’ terms discussed directly below. Their values 
are ad, = 0.165, ag = 0.335, B, = 0.0455, and B = 
0.300. 

The first term of the expression for L,; represents 
the quasi-steady lift calculable from steady-state air- 
It is the lift that would exist were there 
no effects induced by the wake of the airfoil. The last 
two terms represent the effects of the wake expressed 
There is created in the 


foil thee ry. 


in terms of ‘‘lift deficiencies.”’ 
wake, associated with each increment of quasi-steady 
lift, a vortex system whose net effect is to subtract from 
the increment a fractional portion, that portion being 
dependent on the distance of the airfoil downstream 
irom the point at which the increment of quasi-steady 
The subtractive term, the ‘‘lift defi- 
ciency,” is equal to one-half of the increment at the 


lift occurred. 


instant of its creation and decreases thereafter to zero 
as an asymptote according to a lift deficiency function 
that is dependent on 7. The second term of the ex- 
pression for L; is associated with the initial growth of 
quasi-steady lift, while the third term relates to all in- 
crements after the starting instant. 

The deficiency function, defined in reference 2, is 
approximated here by a double negative exponential 
aye r/o bg Ber /b 


Fig. 2 gives a comparison be- 


BY 


ANALOGUE COMPUTATION 387 





NEUTRAL POSITION 









ee 


h 
L, 
/ ts j h' ELASTIC 





ab Y 

~ 
‘ie i] / . 
tie o-ox : cb “a | 
a b s 
2—~ 


Fic. 1 The wing section. 

act at the mid-chord and the three-eighth chord posi- 
tions, respectively, and 1/,, a term representing damp- 
The 
values of these quantities, with senses as shown, can be 
3 


ing in pitch, acts in the negative @ direction. 


given as follows:! 


a 
) — 2rb | a0) 4. : | (aye Bir /b + ase" ") 
aT 


“e 
— Bi/b Be/t a a 
arb / [aye' Bi/b) (+ 71) +. dee! 32/b)( +r 7 ] ( 


tween the approximate and exact functions. The 
approximation is consistent with the degree of accu- 
racy involved in other of the aerodynamic assump- 
tions. 

Referring all motion to the elastic center by the 


transfer equation 
h’ =h + b0[(1/2) —a| (5) 


there are obtained, for the components of force L; and 


L2, 
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Lh 10 
L, = 2ab Jg + : + b[(1/2) — a] no — 2b J 40) aa : 
dr dr | ( 


—(Bi/t — (B2/b 
(aye Bi/b)r 4 ve (Bp ry 


and 
] 42 jd?h +o [0/2 d*9 1 dg) & 
2 = 1O- - ) 24) —-@ i (4) 
ldr? dr? drS ; 
L; and MJ, do not change in the transformation. The 
total positive force and moment on the section are 
z = —TL, _ En of Es (S) 
and 
M=- M+ L,b [ (1/2) + a] — L3b X 


[((1/4) + a] + Leda (9) 


The Lagrange equation is now applied to a three- 
Land M 


dimensional wing of constant chord, with 
taken as the generalized forces. 


Equations of Motion 


The coordinates ) and @ are functions of both 7 and 
the spanwise coordinate x // (x is the spanwise distance 
of a wing section from the centerline of the aircraft and 


/is the wing semispan). Thus, 
h = h(r)y,(x/1) (10) 
6 = O(r)¥9(x/1) (11) 


y, and y, characterize the fundamental mode shapes in 
bending and torsion normalized to unity at the wing 
tip. Asis the practice in flutter work, the mode shapes 
are assumed to remain as they are at zero velocity 
throughout the flight range. (7) and @(7) are taken as 
the Lagrangian generalized 
venience they are referred to hereafter as / and 6, with- 


coordinates. For con- 


x i CAL 
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lh dé 
(0) + d[(1/2) — a] — (0)¢ X 
Ir dr f 


bt a 
. —(Bi/b) ( i) Bo/b . 
ae arb ff (aye ( b) (r T he de 32/b r ) Xx 
0 


j d6 ad*h +ol(1/2 d*6 | 
ldr 2 7U\ 7 dr,2' 


dr," “a 


dr, 


and characteristic of physical properties mass, max 


static moment, and moment of inertia 


ey 
i ¥n7p'd(x/1) 
0 


m = 
a 
o- / ynYop'(€ — a)b d(x/1) I 
0 
a 
[= / ve7p'r?b? d(x/1) 7 
0 


In the last three integrals, p’ is the spanwise line density 


of the wing. In the last integral, r is the radius of gyr, 


tion of a wing section about its elastic axis. 


By consideration of a wing section and integratio: 
along the span, the kinetic and potential energies oj 


the semiwing may be written. Thus, 


r l } (s) ‘ l IT (“) 415 (*)(2 ) 
=> m u \ 
2 dr 2 dr dr af 


V = (1/2)1 moy7h? + (1/2)1 [e0978? (19 


The generalized force for displacement in the h coor 


-" 
—; = l / Lyp d(x /) (20) 
ro 


and for displacement in the @ coordinate is 


ny 
if My, d (x/l) (21 
0 


The equa 


dinate is * ® 


G2 = 


Lagrange’s equation may now be applied. 
tions of motion become 


0 aoc tie f their time dependence ah do “I \ 
out indication of their time dependence. m— + S— + muah = J Ly, a( 99 
The following spanwise integrals are defined: char- dr? dr? J0 l 
acteristic of mode shape ; 
icteristic O I 40 @h ; , - , 
>] + § - Iw, = My, d | - “) 
is — | y,7d (x /) (12) dr* dr? J 0 l 
0 
— " Substituting for L and \/ and considering a wing ol 
A2 = . yo'd (x/1) sani uniform chord, Eqs. (22) and (23) are expanded and 
‘a written in a form convenient for transfer to the analogue 
= " Ynyod (x/1) (14) computer. Thus, the expansion of Eq. (22) yields 
deh dh d°6 ; d6 ad 
(m + 1b?) = (—2rb\.) + h (— mw,”) + 72 (rab*\2 — S) + F [—2rb2(1 — a)A3] + 
dr? at ar° aT 


h 
(0) + [Qi 
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er 
oni b) ) (Bo/b 
nb f [aye “PV? “ — ™) + ase : 
- 9 


d 
6 (—2bdz3) + 2arb (4a.00) + Mi 
cf 


dé br ))T 
YD) — a| A3 ] (0) flare Bi/b + ave B2/b 1) ‘. 
at 


d°6 | 
+ d[(1/2) — a]aAsz of 
dr’ 


dé dh 
+A 


dr ¢ T° 


dt} (24 
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SOLUTION OF 


ind the expansion of Eq. (23) yields 
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"0 { ' ) 2 = fO7h3 >) — ’ f | Pe Irb? a ' d*h 
— | + nb*[(1 8) + a?] A. j= ,<70'A [1 2) a}A2} ~ et — We + 2rb?[(1 eo) a| X24 dr? 
— - dh d9 
rab?hy — S) + {2mb*[(1/2) + a] rs} — 2ab? [(1/2) + a] (4. (0) + As (0) + B[(1/2) — a] We ~ (Oy X 
dr { dr dr 
[aye (Bi/b)r 4 Ave ( Be 1) —2rb?[(1 2) 4. al / [aye (Bi/b) (7 1 (oe 3./b) (r x 
J0 
ee ee. 
5 Ag A; » II 2) — Alhe * aT (2)) 
{ dry ’ dr? ’ dr,24 stile ; 


In the two-dimensional case, the quantities m, S, 
and J take on the meanings of the mass, mass static 
moment, and moment of inertia, respectively, of the 
unit section under consideration. )j, As, and A; become 
unity. In the case of the three-dimensional tapered 
wing, the exponential terms in the lift and moment 
that occur under the time integrals cannot be removed 
irom the spanwise integrals involved in the generalized 
forces. The analogue approach to this complication is 
discussed in Section IV. In Section II, the electronic 
analogue solution of Eqs. (24) and (25) is presented. 

BY ANALOGUE COMPUTATION 


SecTION II—SOLUTION 


The word ‘“‘analogue,’’ when used with respect to 
physical systems, most probably brings to the mind 
of the reader certain component analogies in the vari- 
ous branches of physics. Thus, in some applications, 
the mechanical spring and electrical capacitor, the 
viscous damper and electrical resistor, and the rigid 
body having mass and electrical inductor are considered 
as three pairs of analogous components. They are 
shown to be analogus by recourse to the common mathe- 
matical description of the mechanical and electrical sys- 
tems that they form. This direct component-to-com- 
ponent approach to analogue computation has certain 
advantages, which will be mentioned in Section IV of 
this work. McCann and his associates® and others 
have made successful use of the method. However, 
there are serious difficulties in that a basically new prob- 
lem requires a major application of effort and ingenuity 
to invent a suitable analogue, with no assurance that 
such an analogue can be devised at all. 

The main stream of analogue computation philosophy 
has taken a somewhat broader turn. Rather than 
seeking analogies between components or groups of 
components, it has been attempted to construct units 
that will perform directly, in some physical manner, 
the mathematical operations involved in the mathe 
matical descriptions of the systems to be reproduced. 
Thus the integrator, the adder, and the multiplier have 
become fundamental terms in the computer nomencla- 
ture. It must be mentioned at this point that the dis 
tinction drawn between component and operational 
analogues is not sharp. Thought will show that there 
are situations in which the two are indistinguishable, as, 
for example, in the representation of certain servomech 


anisms in which the original physical components 
carry out such operations as integration. 
portant, the most efficient analogue means often stem 
from a combination of the two concepts. 

Returning to the operational philosophy of compu- 
tation, the first mechanical device for performing the 
basic task of integration was described by Kelvin’ in 
1876. 
good many of the essential ideas of the operational form 
The device described was a 


More im- 


Indeed, in his paper of that year, he set down a 


of analogue computation. 
ball and disc integrator that caused a certain shaft rota- 
tion to be equal to the integral of a variable and con- 
trollable distance in the device with respect to the rota- 
tion of a second shaft. The well-known M.1.T. Differ- 
ential Analyzers* use a more accurate version of this 
unit as their means of integration. A great stride for- 
ward was made in recent years when feedback amplifier 
and servomechanism design reached a point where 
accurate electronic integration, addition, and multi- 
plication became feasible. The outgrowth of this ad- 
vance was a inost desirable product, an accurate and 
flexible analogue computer of reasonable size and cost. 
The REAC, used in the present work, is one example of 
such a computer and is now described in detail. Its 
description will serve to present many of the techniques 
of computation. 

Integration with respect to time, addition, and in 
version (change in sign) on the REAC are accomplished 
by d.c. feedback amplifiers. Fig. 3 is a photograph of a 


typical amplifier. The operations performed by these 





Fic. 3. Integrating, summing, and inverting amplifier. 
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Fic. 4. 


three units are shown in Fig. 4. The inputs indicated 


are applied voltages. Voltages represent the variables 
here as compared to shaft rotations in the M.I.T. Dif 
ferential Analyzers. The value of the voltage supplied 
as an initial condition to one of the integrators of the 
circuit to represent an initial value of a dependent vari 
able determines the scale factor of the problem. Pas 
sage through an amplifier reverses sign. There are seven 
input connections to each amplifier, two inputs having a 
multiplication factor of 10, one having a factor of 4, and 
the remainder having a factor of 1. Fig. 4 also shows 
the schematic representation of the action of a poten- 
tiometer. The basic REAC assembly consists of seven 
integrators, seven adders, six inverters, and 30 potenti- 
ometers with associated power supply, patching equip- 
ment, and Brush-type recorder. Multiplication, divi- 
sion, resolution of vectorial quantities along coordinate 
axes, and trigonometric operations are carried out by 
high-speed servo units. There are also available in- 
put-output tables (for the introduction of special func- 
tions and the plotting of one dependent variable against 
another) and other special equipment. Fig. 5 is a pho- 
tograph of a REAC installation consisting of a basic 
assembly with a six-channel recorder (on the left) and 
a servo cabinet with power supply (the two cabinets on 
Viewed as electronic equipment, the ac- 
Potentiometer and 


the right). 
curacy of the components is high. 
other resistors are held to within 0.1 per cent of design 
value; integrating condensers, to within 0.25 per cent. 
The amplifiers are flat within '/. db. from 0 to 1,000 


cycles per sec. The significance of these figures with 


SCIENCES—JULY, 1950 
respect to the accuracy of solution is discussed jn S¢. 
tion ITT. 

To establish technique, consider the simple problep 
given by the differential equation 


(d*y/dx*) + 20.00 (dy/dx) + 4.14y = 0 — (9 


with the initial conditions y = 0, dy/dx = 0.5 whey 
x =0. Letting REAC time be x, the independent vari. 
able, Fig. 6 gives the analogue setup. One means oj 
tracing the circuit is to isolate the highest order deriya 
tive 

[(d*y/dx?) = —20.00 (dy/dx) — 4.14y| (7 


to assume the existence of —dy/dx at the output of th 
left-hand integrator and to work to the point where that 
assumption is justified by the fact that the input to the 
integrator is d*y/dx* as given by Eq. (27). Applica 
tion of the initial conditions to this circuit will give the 
response of the analogue which can be translated 
directly into the response of the original system. 





Fic. 5. REAC installation 
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Fic. 6. Solution of Eq. (26). 
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FIG. 7 


RC filter and integrator equivalent. 


The ramifications of this technique are both broad 
and extremely interesting from many points of view. 
However, the object of this work is to investigate a spe- 
cific application, and discussion is limited to points 
entering the problem at hand. 

Turning now to the aeroelastic system, one aerody- 
namic-electronic equivalence must be demonstrated 
before the analogue scheme can be devised. Consider 
the simple resistor condenser network and its REAC 
counterpart shown in Fig. 7. The differential equa- 
tion describing both is 


RC(dE» dt) + Eo = E, (28) 


This is easily derived for the RC network and may be 
verified in the integrator case by recognizing that the 
input to the integrator is — (/4)/RC) + (£,/RC) and that 
this quantity, after integration, is —/». Thus, re- 
calling, the sign change through the integrator, 


a | 
-K = —- / [(£,/RC) — (Bo/RC)| dt = (29) 


which, when differentiated and rearranged, is Eq. 
(28). The action of these networks may be found from 
a direct solution of Eq. (28) but, for the present pur- 
pose, is most lucidly seen by means of the weighting 
function concept.’’* The weighting function, W(t), is 
defined for a linear network whose components do not 
vary with time. It is the response of the network to a 
unit impulse input. The weighting function is simply 
related to the response of the network to a unit step 
input, being the derivative of this response. The 
Weighting function of the network under consideration 
is (1/RC)e~’®©, as may be seen by taking the Laplace 

* This concept is, essentially, an expansion on the theory of 
Duhamel’s integral as applied to electronic circuits. 
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transform of Eq. (28) with /, equal to a unit impulse 


and rearranging. This leads to the equation 


Ey = (1/RC) }1/[p + (1/RO)] } (30) 


The stated result is obtained by inverse transforma- 
tion. Any input may be conceived of as a series of 
impulses of proper magnitude imposed successively at 
extremely small time intervals. The linearity of the 
network allows the output to be expressed as the sum 
of the outputs in response to the individual impulses. 
Thus, referring to Fig. 8, the total response of the sys- 
tem at a time ¢ is the summation of the responses to 
the impulses at f;, as 4; goes from 0 to ¢, weighted accord- 
ing to the decrease occurring in the time (¢ — ¢,).. In 
general, stipulating a zero output at /¢ = 0, 


~~ 
=f 
and in this case 
~~ 
Fo(t) = / (1 RC)e™ “ 
0 


A more rigorous treatment of this point will be found 


F(t) Wt —_ tL), (dt, 


W)/RC FB (t)dt; (32) 


in reference 9. 

The mathematical similarity between Eq. (32) and 
the last ‘‘lift deficiency’ terms of Eqs. (24) and (25) is 
The physical analogy may also be drawn if 
The input function 


apparent. 
a convenient construct is used. 


equivalent to “,—namely, 





























dé 4 d*h + alc /2 in a6 
An ~ Am ) 2)-> ” 
dt dr," . dr; 
E,(t,) 
—_ 
*y | 
t, t 
INPUT 
E 
° E(t) 
RC 
t, t 
OUTPUT 
Fic. 8. Weighting function formulation 
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Fic. 9. REAC scheme for lift deficiency terms. 

may be considered an impulse of ‘“‘lift 
deficiency’’ whose effect on the wing at coordinate 
7 is dependent, exponentially, on the distance from the 
wing or coordinate (r — 71). The effect of the defi- 
ciency associated with the initial growth of quasi- 
steady lift can be reproduced by an appropriate initial 


times d7, 


voltage at the network output. 
To summarize the analogy formally, the solution to 


Eq. (28) is 


> " l ar 
E(t) = Eo(O)e~ "8° + / (_-) ge ¢ - 8 Rd 
J 0 


(33) 
This may be verified by substitution. Setting 
i=r, &=s7, I1/REC = B; b 
E)(0) = 2rba; X 


1,00) + Aw ©) + b((1/2 1, 2 @! 
A,G\U) -- Am ;— @ 2)—-@ An (QO) 
a a -"s 


l 
( ) E,(t;) = 2rba; X 
R¢ 


(, dé Ly dh + at(L/2 
An : Am - ) 2) 
\ dr dr; 


d*6 | 


— aja, 
| dr,2) 


it can be seen that the last two terms of, say, Eq. (24) 
are represented by the setup given in Fig. 9. 

It may be noted that if this analogy is carried out 
with the RC network, it can be classed as an example 
of the direct component(s) to component(s) philosophy 
of computation, while, if carried out by the integrator, 
it can be considered a mathematical operation analogue. 

Returning to the complete analogue scheme for 
Eqs. (24) and (25), directly below are given two sets 
of data for two cases of flutter (one two-dimensional, 
one three-dimensional) with constant chord. Case | 
is taken from reference 10 and represents a relatively 
slow speed transport airplane. Case 2 represents a 
high-speed fighter. 


Case 1 
5b = 1,6, = 6it. 
a = —0.40 


AERONAUTICAL 


SCIENCES—JULY, 1950 


mt = 12.59 
os = Zoe 

J = 3.14 
w, = 540/U 
On, = 135 U 

Case 2 

p> = £6 = 3:15. K. 
a = —0.36 
m = 152.4 
5S = 300 

J = 600 


we = 100.2/U 


A, = 0.201 
he = 0.356 
3 = 0.185 


The scheme for the first case is a simplification o 
that of the second and is not presented. The simpli 
fication is brought about by the fact that the ) values 
become unity and the inputs equivalent to /, becon 
equal in both equations of motion. The scheme in th 
second case is given in Fig. 10 for an independent vari 
able scale such that | sec. of REAC time is equal to 10) 
units of r. The potentiometer settings are listed it 
Table 1. Potentiometers | and 
forward velocity. Fig. 10 may be traced through in 
manner similar to Fig. 6 by assuming the existence oj 
d*h/dr* and d*6/dr* at the outputs of summing ampli 
fiers 9 and 10 and working to a point where this as 
The next section will present th 


2 are used to van 


sumption is justified. 
results of machine solutions for the two flutter cases wit! 


numerical comparisons. 


TABLE 1 
Potentiometer Settings 


Potentiometer * Potentiometer * 


1 variable 12 0.799 

2 variable 13 0.214 

3. 60.0216 14 0.482 

4 0.0547 15 0.455 

5 0.106 16 0.505 

6 0.0175 7 0.750 

7 0.0736 IS 0.455 

8 0.274 19 0.510 

9 0.107 20 0.756 

10 0.825 21 0.250 

11 0.262 22 0.115 
* Nore: These readings are uncorrected for potentiometer 
loadings. They represent the values that will satisfy the equ 


tions of motion and differ slightly from final machine settings 


SECTION III—-PRESENTATION OF RESULTS 


Figs. 11 and 12 give a series of responses of the tw 
analogue aeroelastic systems, runs being taken at vat 
ous forward velocities with initial disturbances in / aml 
6. Fig. 11 relates to Case 1; Fig. 12, to Case 2. The 
velocities are closely spaced about the flutter point, tht 
presentation in the first case being restricted to thal 


region. 
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Fic. 10 


Consider the question of flutter velocity and fre- 
quency. Reference 10 gives, for Case 1, calculated 
values of flutter speed and frequency of 832 ft. per 
sec. and 56.5 rad. per sec. Carrying through the 
process of reducing and dimensionalizing the results 
shown in Fig. 11, there are obtained values of S22 ft. per 
sec. and 55.9 rad. per sec. Fig. 13 gives the results of a 
numerical flutter analysis carried out for Case 2 plotted 
in the form of roots of real and imaginary equations 
The flutter speed and fre- 
quency predicted are 580 m.p.h. and 23.4 rad. per sec. 
The 
accuracy achieved with respect to the flutter point is, 


against forward velocity. 
Fig. 12 yields 582 m.p.h. and 23.4 rad. per sec. 


then, satisfactory, being all that can sensibly be re- 
quired in the light of the inaccuracies introduced in the 
formulation of the problem and the use to which the 
results are put. 

The 


pen traces give the total response of the wing in / and 


Consider now the situation at any velocity. 


# (and functions of these variables if desired) to a given 
disturbance. The characteristics of the total trace as 
gathered by inspection are sufficient to give a working 
understanding of the aeroelastic properties of a single 
system. If deeper analysis is desired, it is possible to 
evaluate the modal properties of the response. There 
is, in general, one oscillatory mode that persists in al- 


most pure form after some time. At the higher air 





Analogue for Case 2 


AMP9(5) 





» 


speeds, this is likely to be the mode that finally becomes 
The 
teristics of this mode are directly available from the 
Fig. 14 plots the fre 


negatively damped and causes flutter. charac 


later portions of the records. 
quency and damping of this major mode in Case 2 
The characteristic amplitude ratio 
The 
latter is best taken from a run using higher chart speed 
Means for the study of the 


against velocity. 
and phasing between variables can also be read. 


than that used in Case 2. 
other modes exist in the form of the variable initial 
conditions of the problem which can be used to suppress 
the dominant mode during the initial portions of the 
runs and emphasize a mode of interest. 

Before concluding this section, it is in order to com 
ment briefly on the question of the accuracy of the 
REAC in general and with respect to the present prob 
lem. <A broad discussion of the subject of the errors 
of machine computation has been given by von Neu 
mann and Goldstine.'' The concept of the “mathe 
matical stability’ of a system is introduced by these 
authors, it being a property indicative of the sensitivity 
of a system to small errors in formulation of its mathe 
matical description and the carrying out of the proc 
esses stipulated by the description. Considering only 
the latter operation, a set of equations exhibiting a high 
show “con 


“mathematical stability’’ will 


that is, it will not tend to describe 


degree of 
tinuity’’ of solution 
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U=550mph 














U= 575mph 








U=562mph 
h 








U=61Omph 








Sample runs, Case 2. 


Pic. 12. 


AERONAUTICAL 


SCIENCES—JULY, 1950 

a physical system far diffferent from that intended in the 
face of small variations in parameters and errors in 
mathematical operations. It is clear that the ability oj 
an analogue computer to solve a physical problem ac. 
curately depends on both the “‘mathematical stability 
of the system involved and on the accuracy of the opera 
tions carried out by the computer. Most engineered 
systems display a high degree of ‘“‘continuity”’ of soly. 
tion, since they are designed to operate properly jy 
spite of variable characteristics of components, Th 
REAC has shown consistently good results in the ap. 
alysis of such systems. 
tem, the “‘mathematical stability” can be investigated 


In the case of an arbitrary sys. 


in a number of ways, among them being mathematical 
error analysis of the equations (bringing with it, un- 
fortunately, equations more difficult to solve thay 
those originally considered), machine error analysis 
(carried out by introducing known variations in com- 
ponents and noting deviations in results), and the 
means tacitly employed here of inferring ‘‘stability 
from a number of numerically checked sample problems 
run with machine components chosen at random 
There is no indication in the results of this work that the 
nature of aeroelastic systems makes them markedly 
lacking in ‘“‘mathematical stability.” 

.The next section will discuss the extension of the 
methods described to more complex aeroelastic prob 
lems. 


SECTION IV—APPLICATION TO MORE COMPLEX SYSTEMS 


Some consideration has been given to the extension 
of the present work to more complex aeroelastic sys- 
tems. 
ample, requires special treatment. 
because the aerodynamic wake effects are functions 
of chord lengths of translation in the flight direction, 
a quantity that is a variable along a tapered wing for 
Mathematically, as was 


The tapered three-dimensional wing, for ex- 
This need arises 


a given aircraft translation. 
mentioned, the difficulty is encountered in attempting 
to remove the “‘lift deficiency’? nondimensional time 
integrals from the spanwise integrals given by Eqs 
(20) and (21). As is the practice in numerical flutter 
analysis, the initial approach has been to divide the 
wing into a number of constant chord sections (this for 
aerodynamic purposes only) and to carry out the spat 
wise integration by summing the integrations ove! 
constant chord sections. On the REAC this can be 
accomplished by an array of equivalent RC filters (Fig 
9). Since such use of integrators is wasteful, in the 
tapered wing problem it is advantageous to use a grou 
of true RC networks. This allows many constant chord 
sections to be taken without an unreasonable amout! 
of equipment being required. 

The question of structural damping can be handled 
on the REAC in almost any way the analyst desires 
Viscous damping can be included, of course, by trivial 
modification of Eqs. (24) and (25) and no modification 
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SOLUTION OF 


of equipment. Damping always in phase with ve- 
and proportional to amplitude of motion can be 


locity ‘ 
d by the use of a few simple relays. 


produce 
The consideration of additional degrees of freedom 
s two additional integrators and related equip- 


involve 
There seems to be no 


ment per degree of freedom. 
intrinsic reason why accuracy should deteriorate with 
these additions, this again being determined by the 
change of “mathematical stability’’ of the system as a 
function of the additions. 


As has been mentioned, McCann and his coworkers 
have solved various problems with the use of passive 
network, nonoperational analogue techniques. Among 
these problems is that of beam vibration in torsion and 
bending. The analogue is carried out by breaking the 
beam into lumped masses and inertias connected by 
massless springs and representing the lumped elements 
by suitable passive networks. This is a difference 
equation technique and can be extended to solve the 
partial differential equations of aeroelasticity. Mc- 
Cann has recently made known work in this di- 
rection for the problem of gust entry. It is possible to 
attack the partial differential equations with REAC 
techniques by considering the lumped system and asso- 
ciating a new time dependent variable with each 
element. 

It is too early to draw a comprehensive picture of the 
techniques and equipment that will prove most suitable 
for the various aeroelastic problems, though one might 
guess that an operational computer supplemented by 
proper passive networks would form a most versatile 
unit. Clearly, further work is needed in the direction 
of rigorously exploring the accuracy and capabilities of 
the analogue computing techniques. 
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othesis for Skewness of the Probability 


Density of the Lateral Velocity Fluctuations 
in Turbulent Shear Flow’ 


STANLEY CORRSINt 
The Johns Hopkins Unwersity 


SUMMARY 


It is proposed that the observed skewness in the probability 
density of transverse velocity fluctuations in a turbulent shear 
flow is due to the corresponding lateral gradient in lateral turbu- 


lence intensity. A simple estimate is deduced for small skew 
ness in terms of the turbulence level, its gradient, and a char- 
acteristic length, A. An order of magnitude comparison is made 
with experiment by computing A and comparing it with the 


dimension of the shear region. 


NOTATION 


y = lateral coordinate (perpendicular to mean flow and to 
line source) 


l = time 

v = lateral velocity fluctuation 

0 = mean velocity 

w = 0/0 

6(w) = probability density of w(t) 

w’ = Vw 

wy’ = turbulence level at the heat source location, y = 0 


U HAS BEEN FOUND EXPERIMENTALLY that the trans- 
verse temperature distributions close behind fine, 
straight, heated wires set perpendicular to both mean 
flow and direction of largest gradients in two turbulent 
shear flows are skew.'~* The thermal wake close behind 
a line heat source in isotropic turbulence is symmetri- 
cak.* ° Taylor's theory of ‘diffusion by con- 
tinuous movements,’’® it is evident that in the latter 
flow this transverse temperature distribution is simply 
the probability density of the lateral velocity fluctua- 
tion, v(t). For shear flow this correspondence should 
still be at least approximately true. 
is an attempt to explain this apparent skewness in the 
probability density 6(w) as being a result of the lateral 
lateral Dryden has 
that this measured skewness is much 


From 


The present note 


gradient in turbulence level. 
pointed out? 
greater than that which could be predicted from gra- 
dients in mean velocity in the wake or from diffusivity 
gradients due to differences in turbulence level within 
the wake. 

The probability density must satisfy the basic rela- 
tions, 
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* I should like to acknowledge the helpful criticism of Dr. 
Leslie S. G. Kovasznay and to thank Dr. G. B. Schubauer for 
making available the unpublished data of the measurements re 
ported in reference 1. 
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f- 6 dw = l (] 
f- wb dw = 0 (2 
f- wd dw = ay? (3 


The ‘‘skewness factor”’ is defined as 


‘ Wo’ l ile 
S= =; w*§ dw (4 
(W 9° Wo « x 


(Ow’ /Oy 
and a characteristic length A can be deduced with the 


A simple estimate of S in terms of wy’, 


restrictions of small turbulence level gradient and small 
We define n 


*0 
i. = / 6de: n. = / 6 dw 5 
F : 70 


n_+n,= 1 6 


skewness. and 2, by 


Thus, 


and, for small skewness, each should differ only slightly 
from the value '/.. We characterize the negative and 
the positive velocity fluctuations by the quantities w 


and w,, respectively : 


-_ ‘ 
tw? = / wh dw; niws,? = | wh dw (7 
, 0 


Thus, 
Pt niws? = wo’ (o 


and w_ and w, should differ slightly from wo’. Actu 
should be turbulence levels character- 


at some lateral dis- 


ally, w_ and w, 


istic of the turbulence intensities 
tances from the line source (y = 0). Since negatively 
directed velocities originate in the y > 0 region, we ex- 
pect w_ to be a measure of the turbulence level at y = 
+A (say); conversely, w, should be the turbulence level 
at some distance y = —A. The positive and negative 
lengths are taken of equal magnitude, since this 1s 4 
small perturbation calculation. In effect, we are using 
the first part of a series expansion for the turbulence 
0, the location at which 


Hence, 


level distribution, about y = 
we are interested in the skewness of 6(w). 
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VELOCITY FULUCTVUATIONS 
wW = — wy’ (Ow’ Oy)oA t 


’ , (9) 
Oo, = wo — (Ow Oy)oA f 


These may, 10 effect, be considered as the definition of 


a characteristic length A. 

In order to get a quantitative estimate of S, some 
typical form of 6(#) must be assumed. A small devia- 
tion from Gaussian can be constructed by recalling that, 
for zero Skewness, 

ae 1 w | 
= exp ) — * 


/ « y= 
Vv 27 wo z Wo 


Hence, for the present estimate, we choose 


fr -m < we —- C¢ 
| j (w + c)? ) 
b(w) = : exp< — , -> 
V/ 2t wo (V/2@ — e)?f | 
, > (10) 
for -C< wc @;: 
: | j (w + c)? } 
6(w) = -exp< — : -> 
\ Ir wo { Vv r.4 Wo +- { 


is illustrated in Fig. 1. The restriction to small skew- 


ness requires ¢ wo’ << 1 and €/wo’ << 1. The quan- 
tities cand e are defined such that c, e > 0 corresponds to 
S< Qand vice versa. 

Insertion of Eqs. (10) into Eq. (2) gives the relation 


between c and e: 


e = (1/2) Y 1c (11) 
From Eqs. (5), 
n_ = 4 l + 2=— *) c| 
2% V 2n wo’ \ r { 
(12) 
=-J1- (2-Z)cl 
Z { NV 2r wo’ 2 f 


which obviously satisfy Eq. (6). 
From Eqs. (7), neglecting terms in c? and higher pow- 


ers, 

, - 

Wo , ” 9 , 

i w.* = ( — ) V 2 wo C 
2 T Ss 

(13 

) 

Wo l ) 5 , 

Ny Ww 3 ( = ) V2r wc } 

2 ve 8 


Division of Eqs. (13) by the corresponding members 
ot Eqs. (12) yields approximate expressions as fol- 
lows: 


(14) 


Comparison of these with Eqs. (9) shows that 
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Finally, with c and ¢e evaluated, an approximation 


for skewness results from substitution of Eqs. (10) into 
Eq. (4), after higher order terms in c and ¢ have been 
neglected, 


Ss = —A(A wo’) (Ow’ OY)o (16) 


Rough comparisons with experiment are afforded by 


the measurements of references 1 and 3. Since A is 
relatively vaguely defined and since all of the other 
quantities in Eq. (16) have been measured, the check 


will consist of solving for A in Eq. (16) and seeing 
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whether its order of magnitude is consistent with the 
size of the turbulent shear flow region: 

(1) The line source was set parallel to the wall in a 
turbulent boundary layer (Fig. 2). The experimental 
results were: 


Ow’ 


wo’ = 0.028; ( ) = —0.0070; S = 0.38 
0 


bf 
Then, Eq. (16) gives 


A = 0.38 cm. 


IVEY, 13950 


A = 1.1 em. 


Both of these lengths appear to be of reasonable order 
of magnitude, in spite of the rather divergent properties 
of the two flows. 
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(II) The line source was set perpendicular to a 
radius in a turbulent jet (Fig. 3). The experimental 
results were: 

Ow’ 


) = 0.05; S = —0.10 
Ov /o 


Eq. (16) gives 
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Longitudinal Stability of Autopilot- 


Controlled Aircraft 


ANDREW VAZSONYI* 
Naval Ordnance Test Station 


ABSTRACT 


The design analysis of autopilots for subsonic and supersonic 


aircraft is presented 
in a concise operational form 


The aerodynamic equations are written 


with the aid of a new set of aero- 


dynamic coefficients which has immediate bearing on the dynam- 


ics of the aircraft. 


Study of the uncontrolled aircraft is pre- 


sented first, then the analysis of a proportional type control is 


given. 


The final analysis contains a properly tuned computer 


between the sensing element and the autopilot, which is assumed 


to have a simple time lag.s Numerical examples are presented for 


both a subsonic and a supersonic aircraft. 


The analysis is based 


on the transient behavior of the system and is further sub- 


stantiated by transients obtained with the aid of an analog com- 


puter 


NOTATION 


40, Ai, A = 
Bo, By, Bo, B = 
c = 
1 

m = 
Mu, Me, M;, My = 
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coefficients of cubic equation 

coefficients of fourth-order equa- 
tion 

capacitance of condenser 

gravitational constant 


moment of inertia of aircraft i 
pitch 

mass of aircraft 

aerodynamic coefficient, pitching 
moments per unit of u, w, w, and 
q 

aerodynamic coefficient giving 
pitching moment per unit of 
change of elevator deflection, 6 

differential operator 

angular pitching velocity 

velocity of aircraft along longi- 
tudinal axis 

change in velocity along longi- 
tudinal axis 

velocity of aircraft along normal 
axis 

change in velocity along normal 
axis 

resistance 

time 

time constant in Eq. (48) 

time lag of autopilot 

anticipation time constant in Eq. 
(94) 

voltage in Eq. (120) 

aerodynamic coefficient, force 
along longitudinal axis per unit 
of u, w, and g 

aerodynamic coefficient, force 

along normal axis per unit of 

u, Ww, and q 
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4 


ay, ae 


wy 


aerodynamic coefficient, force 
along normal axis per unit of 
change of elevator deflection, 6 

change in angle of attack of air 
craft measured from steady 
state angle of attack 

dimensionless coefficients for cubic 
charts |Eqs. (132) and (133)] 

change in elevator angle measured 
from steady-state position 

control damping coefficient (di 
mensionless ) 

aerodynamic coefficient or dimen- 
sionless wind-tunnel damping 
ratio [Eq. (11)] 

dimensionless damping ratio of 
fast oscillations 

dimensionless damping ratio of 
slow oscillations 

aerodynamic coefficient or dimen- 
sionless ‘‘gyroscopic’”’ damping 
ratio 

aerodynamic coefficient or dimen 
sionless phugoid damping ratio 

change in attitude of aircraft 
measured from steady-state atti 
tude 

steady-state attitude 

aerodynamic coefficient, change 
in angle of attack per change of 
elevator angle |Eq. (27)] 

dimensionless follow-up ratios 

dimensionless aerodynamic co 
efficient in Eq. (19) 

voltage from position gyro 

aerodynamic coefficient, time lags 
defined by Eqs. (24) and (63) 

undamped natural frequency of 
fast oscillations 

undamped natural frequency of 
slow oscillations 

dimensionless coefficient on cubic 
chart 

undamped natural “gyroscopic” 
frequency 

undamped natural wind-tunnel 
frequency 

aerodynainic coefficient, fre- 
quency defined by Eq. (57) 

undamped phugoid frequency 


INTRODUCTION 


remy FOR CONVENTIONAL AIRPLANES are de- 
veloped and designed without making elaborate 


dynamical stability studies. 
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Recent experience with 
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46 | COMPUTER o AUTOPILOT | re 
o/8 =Le(P) |e La) | 





INTERIOR DYNAMICS 
§/0=L(plLa(p) =£i (Pp) 





8 | EXTERIOR DYNAMICS 
7 L 8/6=£ (p) 


Fic. 1 Notation and block diagram. 





pilotless aircraft and, in particular, guided missiles indi- 
cates that, for high-speed aircraft, servo systems can- 
not be developed unless a thorough understanding of 
the dynamics of the aircraft and the servo system is 
available. 
a dynamical study for both a conventional and a super- 


The purpose of this paper is to present such 


sonic aircraft and also to point out the similarities and 
differences between the two cases. A moving aircraft 
controlled by an autopilot presents a complicated prob- 
lem in dynamics. In order to organize the analysis, it 
will be necessary to examine this problem of dynamics 
in a systematic manner and outline the method of ap- 
proach to be employed. 


(1) Method of Approach 


The motion of the aircraft can be represented by three 
We will restrict 
our discussion in this paper to the so-called longitudinal 
In this case, the lateral linear motion and the 


linear and three angular coordinates. 


motion. 
rolling and yawing motion are assumed to be zero. 
When the motion is restricted to small perturbations, 
the motion can be described with the aid of a system of 
three ordinary linear differential equations connecting 
the perturbations of the forward velocity (#), normal 
velocity (w), angular pitch velocity (q), and the ele- 
vator deflection angle (6). 

Fig. 1 shows a so-called block diagram of the system. 
The exterior aerodynamics is characterized by the longi- 
tudinal equations of motion. Instead of the variable 
q, its integral 6, the change in attitude of the aircraft, is 
used. 
termining the attitude function 6(/) for any given ele- 
vator deflection function 6(¢). 

As far as the interior dynamics is concerned, it will 


The exterior dynamics is characterized by de- 


be assumed that the aircraft is equipped with a direc- 
tion measuring device (like a free gyro) which senses 
6, the change in attitude of the aircraft. The interior 


dynamics of the system determines the elevator deflec- 
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tion function 6(¢) for any given attitude function 4); 
For purposes of analysis and design, it will be necessary 
to break the interior dynamics into two separate con 
ponents: the computer and the autopilot proper, |; 
will be assumed that a signal proportional to the chang 
in attitude @ is fed into the computer; that the com 
puter determines first, second, and higher order 9 
derivatives and integrals; 
from the computer. 
pilot proper, which in turn actuates the elevators. The 


and that a signal o results 
This signal o excites the auto 


change in attitude @ is connected to the signal ¢ by » 
system of differential equations characterizing the com 
puter, while o and 6 are connected by a system of differ 
ential equations characterizing the autopilot. 


In many servo problems—1.e., torque amplifiers and 


follow-up mechanism—the “‘exterior’’ dynamics js 
simple, and servo theory is restricted to dynamics of the 
“interior.” However, for airplanes this is not the case, 
since even the motion of an uncontrolled airplane pre 
For this 


reason, part of the effort in this paper will be to endeavor 


sents a complicated problem in dynamics. 


to obtain a clear physical picture of the exterior dynam 
ics of the plane. Thus, the longitudinal equations oj 
motion will be examined in great detail, and a physical 


feeling for these equations will be developed. 


The longitudinal equations of motion represent three 
ordinary differential equations for the variables uw, w, 
and qg. As it is more convenient to work with the 
change in attitude (@) and change in angle of attack 
(a), the equations are written with these variables. 
The coefficients of these equations are the mass and 
inertia of the plane and the aerodynamic coefficients 
However, it would be difficult, indeed, to infer the 
dynamic properties of the system just by observing the 
magnitudes of these coefficients. One important aspect 
of the theory presented in this paper is the fact that the 
equations of motion can be rewritten in such a form 
that just by observing the coefficients of the equa 
tions the important dynamical features of the aircraft 
can be perceived. The most important quantities in 
dynamics are frequencies, time lags, and dimensionless 
damping ratios. It is possible to introduce for every 
airplane a set of frequencies, time lags, and dimension 
less quantities which completely describes the system 
These quantities can be computed from the aerody 
namic coefficients with the aid of simple algebra, and so 
they can be considered as a new set of aerodynamic co 
efficients. Two facts are to be emphasized about these 
coefficients: First, they can be computed without the 
solving of any algebraic equations of higher order 
Second, these coefficients do not depend at all on the 
We call them aerodynamic 
coefficients, though the mass and inertia of the aif- 
[As an ex- 


autopilot of the system. 


craft enters into determining their value. 
ample of the use of these coefficients, we point out 
Eq. (23), 
dynamics of the airplane. 


which completely describes the exterior 
In Eq. (23), all w's are fre- 
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STABILITY OF 


quencies and all ¢’s are dimensionless damping co- 


efficients. ] 

After a complete understanding of the exterior dy- 
namics is achieved, attention can be focused on the 
analysis of the interior dynamics and then of the com- 
plete loop. The simplest mathematical assumption is 
to take the elevator deflection 6 proportional to the 
attitude 6. This assumption is unrealistic, since it 
ignores autopilot time lags. (In other words, a sudden 
change of attitude will require a sudden change of ele- 
yator deflection which, of course, is not realizable phys- 
ically.) In spite of this idealization, important con- 
clusions can be deduced from the analysis of such a pro- 
In particular, the peculiarities of 


portional system. 
The crux of the 


airplane controllability become clear. 
matter appears to be that, while an increase or decrease 
in lift force is required to maneuver, there are actually 
no direct means to change the lift of the aircraft. The 
motion of the elevator applies a torque to the aircraft; 
because of this torque, the aircraft starts to change its 
attitude. For the subsonic case, this motion is essen- 
tially nonoscillatory, and therefore no particular trouble 
arises. However, in the supersonic case, the variation 
in attitude is oscillatory. Only after these oscillations 
die will the aircraft achieve a new angle of attack and 
will lift be developed to accelerate the plane up or down. 

The thorough understanding of pure proportional 
control will help to synthesize the interior dynamics of 
the aircraft and aid in obtaining satisfactory perform- 
ance. In particular, it will be shown that, for the con- 
ventional aircraft, a simple proportional-type control is 
In the case of the supersonic aircraft, 
First, the fast 
that is, the ele- 


satisfactory. 
the situation is more complicated. 
pitch oscillations must be filtered out 
vators must not be permitted to follow the fast oscilla- 
tions of the aircraft. Furthermore, as there is not 
sufficient aerodynamic damping available in pitch, it 
will also be necessary to move the elevators in a manner 
that introduces artificial pitch damping. It will be 
shown that to accomplish this requires a much faster 
autopilot than the one used for the conventional air- 
craft. 

The next problem arising is the physical realization 
of the interior dynamics arrived at through the syn- 
thesis. The autopilot proper is usually a complicated 
system, and its dynamics can be described only by a 
complicated set of differential equations. As it is not 
the purpose of this paper to go into the detailed design 
study of the autopilot, a drastic simplification will be 
made in this respect. In particular, it will be assumed 
that the dynamics of the autopilot can be expressed by a 
simple exponential time lag. In mathematical terms, 
this means that it will be assumed that the signal and 
the elevator motion are connected by a linear ordinary 
differential equation of the first order. The magnitude 
of this time lag is an indication of the speed of response 
of the autopilot. In the case of the conventional air- 
craft, it can be shown that the autopilot time lag is of no 
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particular importance. However, in the case of the 
supersonic plane, it will be shown that the autopilot 
time lag must be short, or, in other words, the autopilot 
must be fast. 


Part I—SuBsonic AIRCRAFT* 


(2) The New Aerodynamic Coefficients 
The longitudinal behavior of an aircraft is character- 
ized by the following set of partial differential equa- 


tions:! 


m du — X,u — X,w (mW — X,)q + 
dt 
mg cos wf gdt=0 (1) 
dw : P 
—Z,u +m — Z,w — (mU + Z,)q + 
at 
mg sin A fo dt=2Z,6 (2) 
dw dq 
—-M,u — M; — M,w+Tl1 —- Mq=M,6 (3) 
dt dt 


where Eqs. (1) and (2) are force equations along the 
longitudinal and normal axis of the plane and Eq. (3) 
is the pitching moment equation. It will be more 
useful to introduce a, the change in angle of attack of the 
aircraft, and 6, the change in attitude, both measured 
from the steady-state value. These quantities are 
given by 


w/U 
S qdt (5) 


Eqs. (1), (2), and (3) become now, in operational nota- 


a= 


d= 


tion, 
m m 


| (w — xs) pt gost |0 =@ (6) 
m 


—£ ’ —Ze 
4+ Zu u+ l (» + ) an 
m m 
I(v + =) p+ gsin | 6 = 4s 6 (7) 
m m 
—M U 
“u + —(—Myp — Mya + 
" I I 
: —M, ) ; M;. ; 
b? + D100 = 6 (8) 
ees 


Interpretation and numerical values of the above 
aerodynamic coefficients can be found, for instance, in 
reference 1. No discussion of these coefficients will 

* Sections (2) and (3) are applicable to subsonic or supersonic 
aircraft and are included here as an introduction of the subject 
and because these sections are more pertinent to the subsonic 


case. 
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be given here; as for the dynamic stability studies to be 4 ei » P+ 2ewep + we” - 
rsued in this paper, a new set of aerodynamic co- PE 2S HTT 200" o+ “" 
pursued in this paper, a new set of aerodynamic co (mp + 1)(r2p + 1) is sit 
efficients will be used. In order to simplify matters, hw, 25 = 0) (95 ner 
. ia . * — 4) e 
some of the small aerodynamics coefficients of Eqs. the 
(6), (7), aud (8) will be neglected. The justification of Interpretation of this equation will be given in the next than 
this can be obtained a posteriori when an actual nu- Section of the paper. ey 
S ¢ 
merical example is being considered. In particular, the ; : ; : 
. i ; é (3) Interpretation of the New Aerodynamic Coefficients cae 
following assumptions will be made: ts 
' : : (a) If the aircraft is suspended at its center ae 
X,= 2, = My = M,=0, 2=0 (9) gravity in a wind tunnel, obviously the change in ats. | 
) J 5 1 att be ac 
F ; tude of the aircraft equals the change in angle of attag| — 
Let us now introduce the following new aerodynamic ae : 1 § gle of attack in th 
oa . or?=a. Thus, Eq. (22) becomes 
coefficients: the § 
‘ ~ (p? + 26 wWwP + Wy") a + AW "4 = (0) (6, will 1 
wy? = —M,U/I (10) " 
/ Under static conditions, a fixed change in elevator de 
fe — a 9 ame y ° . e ° 
i, = —M,/2 MUI (11) flection gives a fixed change in angle of attack or This 
= — ie 2 ‘ g 
A= M; M,l (12) a = &5 (7 undé 
— il r (Y: : ee scop 
we” = —(Zyg cos 0o/mU) + (Xig sin 6/mU) (13) Hence, d is the ratio of angle of attack change to elevator In 
2tawe = —(X,/m) — (g sin /U) — (Z,W/mU) (14) deflection angle chang . Therefore, d is a static measure be w 
of the elevator effectiveness. 
We = (—LyXw + XuZy)/m? (15) Assuming now that the elevator is locked and 6 = (), 
$ , r . < 96 20, 2c 
2¢ we, = (—Xu — Zy)/m (16) Eq. (26) becomes whe 
. Q. r ir (p* + 2 FuWeP + wy*)a = 0 (28 lags. 
a % = X,,/m (17) give 
. rhis is the equation of a damped harmonic motion so me 
Q, = —X,,/m (18) : : hig type 
k " common in mechanics, circuit theory, and servomecha- A 
> - 1° > r - > > Pos > 4¢ r . ‘ npat: . 
t = mg/UX, (19) nism theory. The moment of inertia of the aircraft wine 
takes the place of the moment of inertia of the pendu- char 
where the w’s and Q’s all have the dimension of fre- lum, while gravity is replaced by the spring-type action scop 
quency (the reciprocal of time) and where {o, A, and £ of the lift forces. Because of the angular motion of the 
are dimensionless. With the above notation, the equa- aircraft, there is an apparent angle of attack between 4 
tions of motion [Eqs. (6), (7), and (8)] become the aircraft and the air stream. This latter angle of ™ 
7 10,/( )] ( ae - attack results in a lift force proportional to the angular 
u = [2)/(p + Q)| (a — Ecos 4:0) (20) . ‘ i , : 
tp + M)I | ’ ’ velocity of the aircraft and yields the damping term of E 
(Pp? + Wawa + we2)a = (p? + 2eewep + wo?)0 (21) the oscillatory motion. From Eq. (28) it follows that 
wy» 18 the undamped ‘‘wind-tunnel”’ frequency of the air 
(p? + 2 wwwP)O + wy?a + Awy?d = O (22) craft and §» is the dimensionless wind-tunnel damping Thi 
‘age ; nis 
: ey a _ ratio of the aircraft. 
Eq. (20) is independent of Eqs. (21) and (22), and it 7 oa? : ; ; ee and 
‘ Se Z : es The above statements give immediate physical 
serves for the determination of « when desired. The . - " ; 
: ange : : meaning to the coefficients w,, (,, and X. 
angle of attack a can be eliminated from Eqs. (21) and ee we ; , 
») + thie Selo salle etatiilie nay tye (b zancester developed the so-called phugoid dil 
22), and the following single stability equation is ob- pe a Pe 
ht? d 5 =e y eq theory of aircrafts’? wherein it is assumed that the 
ained: ne eae ie 
moment of inertia of the plane is negligible and that the 
2 oy 2 aircraft instantaneously adjusts itself a constant 
— (PE + 2Wowop + wo? , urer ult instant neous!) idjusts itself to a constant 
PP + WFuwup + wu a - + angle of attack. Thus, a remains zero and Eq. (21) be 
p- —- 2C Wal - Wa” m n Ww 
2 comes, now, 
Aw,” 6 = O (23) 
; : ; ; (p? + 2owep + we”)? = O (29 
An alternative useful form of the equation can be ob- 
tained by the use of the parameters Again, this is the equation of a simple damped harmonic The 
motion. It can be seen that ws is the (undamped) phu- acte 
(e+ V/¢ ° 1) | (< Ve? 1) goid frequency of the aircraft, while ¢» is the dimension- slov 
= wal Gal ila 5) =OeSe = a la a ; : eh ; 
7 matin . T2 ‘ less phugoid damping ratio. While executing the phv- In 1 
(24)  goid oscillations, the aircraft moves on an oscillating be f 
' ' he di . — witt path. At the peak of the waves, the forward velocity (p 
where 7; and 72. have the dimension of time. ith . : 
: es of the plane is lower than at steady state, and so the 
these parameters, Eq. (23) becomes ‘ ete ace , ou : : 
aircraft's kinetic energy is below average. Because Ol whe 
* Some of these coefficients have been used by other researchers. _ the slowing down, the lift is not sufficient to support the rati 
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STASILITY oF 


ht, and the plane starts to sink. While the plane 


weig 
inking, the potential energy transforms into kinetic 


rgy, and the craft speeds up. At the low points of 
the wavy path, the velocity is high, there is more lift 


is 


ene 


than needed, and the plane starts to rise. Thus, there 
is a periodic exchange of potential energy into kinetic 
energy giving rise to the up-and-down motion. 

Finally, the aircraft is assumed to be flying with 
(This could 


(Cc) 
a constant attitude, and therefore 6 = 0. 
be accomplished, as an example, by carrying a large gyro 
in the aircraft. When a pitching moment is applied, 
the gyro will precess, but the attitude of the aircraft 


willremain unchanged.) Eq. (21) now becomes 


(p? + 26,wap + w,”)a = 0 (30) 


Thisis, again, a harmonic motion; w, is the “‘gyroscopic”’ 
undamped natural frequency, while ¢, is the “gyro- 
scopic” dimensionless damping ratio. 

In most cases ¢, is larger than unity. 


be written, with the aid of Eq. (24), as 


(1 + pr) (1 + pria = 0 (31) 


Eq. (30) can 


where 7; and 72 could be called the “‘gyroscopic’’ time 
lags. The physical meaning of these time constants is 
given by the speed of the nonoscillatory (subsidence) 
type of motion. 

As a summary, let us state that: w, and ¢, are the 
wind-tunnel characteristics; we and ¢» are the phugoid 
characteristics; and w,, ¢, 71, and 72 are the “‘gyro- 
scopic’ characteristics. 


4) Uncontrolled Flight 
0, and Eq. (23) becomes 


p> + 2fewep + “| 6@=0 
Pp? + 26 Wah + wo,” 


In this case, 6 = 


EG + 2 WwP + ww” ( 
(32) 
This equation is a fourth-order differential equation 


and can be rearranged as 


(pt + B3sp? + Bop? + Bip + Bo)@ = 0 = (33) 


where 
Bo = wy*we" (34a) 
By = 2(CuWuwe” + Cowewy”) (34b) 
Bg = wy? + wg? + Af wl WaMu (34c) 
Bz = 26 a@a + SuWw) (34d) 


uncontrolled aircraft can be char- 
one fast and the other 


The motion of an 
acterized by the two frequencies 
slow—and by the two corresponding damping ratios. 
In mathematical terms, this means that Eq. (23) can 
be factored into two quadratic factors. 

(p? + 2e¢,wpp + w/) (p? + WwsP + w;7)0 = 0 (35) 


where w, and ¢, are the fast frequency and damping 
ratio and w, and ¢, are the slow frequency and damp- 
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ing ratio. By comparing coefficients, the following 


equations are obtained: 


Bo = wy*w;” (36a) 
By = 2(Cpws + Fwy) wyws (36b) 
By = wf + wy.” + 45 /f.@.0y (36c) 
Bs = 2(Cyws + Fas) (36d) 


In order to determine wy, w,, ¢, and ¢,, the roots of the 
fourth-order equation [Eq. (33)] must be determined. 
However, for practical application, 


Wy > ws (37a) 
wy? D> 45 Fw yrs (37b) 
Coe > F5as (37c) 


and, therefore, the second and third terms of the right- 

hand side of Eq. (36c) can be neglected. Thus, 

we = V By = wyV1 + (wa?/wu?) + AS wha (Wq/ ww) 
(38a ) 


Furthermore, the second term of the right-hand side of 
Eq. (36d) can be neglected, and, therefore, 


¢, = B;/2w; = B;/2~/B, (38b) 
From Eq. (36a) and Eq. (38a), it follows that 
a, = WV Bo/a; = W Be/B: (38c) 


and, finally, with the aid of Eq. (36b), it follows that 

" | ( B, ‘ ) B, 

a —™ $j = sn 
wy \2w ws 


2/ BoB 
Eqs. (38) are approximate and depend on the validity 
However, for any 


l Bs Bo 
2 Bs \ By 
(38d) 


of the inequalities of Eqs. (37). 
actual numerical example, an a posieriori verification 
is possible. It is believed that for most cases Eqs. (38) 
are sufficiently accurate. 

It will be useful at this point to introduce a numeri- 
cal example. A light fighter plane has the following 
characteristics :* 





0.2 sec. 


One div. = 


Fic. 2. Wind-tunnel transient. 
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X,/m = —0.074; Z,/m = —0.323; M,/I = 
— 0.0364 
(39) 
Xu/m = 0.094; 2Z,/m = —2.28; M,/I = 
—4.15 


All other coefficients are zero. The above figures - 
English units. The velocity of the plane is 252 ft. per 


sec., and its weight is 4,700 Ibs. 


The new aerodynamics coefficients can be easily computed for level flight : 


wind-tunnel characteristics, w, = 3.03 rad. per sec.;  {» = 0.685 

phugoid characteristics, w» = 0.203 rad. per sec.; fo = 0.018 

gyroscopic characteristics, w, = 0.446 rad. per sec.;  (, = 2.68 (40 
7, = 0.442 sec.; Tze = 11.2 sec. 

period of phugoid oscillation = 31 sec. ; 


Thus, the wind-tunnel frequency is comparatively 
high (0.48 cycle per sec.), but it is well damped (Fig. 2). 
The phugoid frequency is much lower (0.032 cycle per 
sec.) and is practically undamped. (Figs. 2 and 3 show 
these transients when a sudden moment is applied to the 
aircraft.) 

Eqs. (38) give the dynamic characteristics of the un- 
controlled plane. Numerically, 


Wr = 4.37 f, = 0.745 | 
w, = 0.141 ft; = 0.141f (41) 


Period of long oscillation = 44.5 sec. 
The assumptions involved in Eqs. (37) should be 
checked for the numerical values considered. 





Fic. 3. Phugoid oscillations. One div. = 4 sec. 





Fic. 4. Transient response of conventional aircraft. One div. = 
+ sec. 


From Eq. (38a), the value of B; can be found as 
Bz = wf = 13.7 (42 
while, according to Eq. (36c), 
By = ws + Ws” + ATF wyws = 13.7 + 0.04 + 
0.26 = 14.0 (43 
The value of B;, according to the first half of Eq. (38b), 
is given by 
Bs = 2wyt, = 6.52 (44 
while, from Eq. (36d), it follows that 
Bs = 2wyfy + 2wsfs = 6.52 + 0.038 = 6.56 (45 
Thus, it can be seen that the discrepancy in B, and B;is 
negligible. (Bo) and B, need not be checked, since there 
is no approximation involved in their determination. 
Fig. 4 shows the transient behavior of the attitude of the 
disturbed airplane when a sudden moment is applied to 
the aircraft. 
(5) Proportional Control 
We will assume that 


6 = po (46 


Eq. (23) becomes 


G = 2 Fw WwP + wn -— saa : ~“)+ dae’ X 


p2 + 26 Hab + wa" 
¢6=0 (4 
or 
(pt + Bsp*? + Bop? + Bip + Bo)é = 0 (48 
where 
Bs = 2(CaWa + Sww) (49 
Bo = wy? + wa? + 4falwWuWa + Amwn” (50 
By = 2W(CyWnwe? + Cowow? + Apwn* awe) (1 


Bo = (Ww?we” + Apw,,"w.") (52 


In order to study the system, the roots of Eq. (4 


should be determined. Equations similar to Eqs. 
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Fic. 5a (top). Fig. 5b (bottom). 


(38) give approximate values of the frequencies and 
damping ratios. Figs. 5a and 5b show the plot of wy, 
®,, fy, and ¢, for the airplane studied in Section (4). It 
will be useful to study these curves carefully. 


When Au increases from zero (B; remains unchanged; 
By, By, and By all increase), wy, ws, and ¢, increase while 
¢; decreases. When Au = 0.46, ¢; = 1, showing that 
at this value of \u the slow oscillations degenerate into 
a slow nonoscillatory motion (called subsidence). 
Thus, it can be seen that, by the proper selection of Au, 
a system can be obtained for which both damping ratios 
are good. An increase in \uv will make the systems 
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tighter, but too large a value of Aw will decrease £,; too 
much. Figs. 6, 7, and 8 show transients of the atti- 
tude (and the derivative of the attitude) of the dis- 
turbed aircraft when a sudden moment is applied. 
The value of \u = 1.0 would probably give the solution 
considered optimum. It can be seen that the aircraft 
will take up a new attitude in about | sec., and there- 
fore it can Le concluded that no extreme speed of re- 
sponse will be required from the autopilot. 























3 
te 
dan - 
Fic. 6. Transient response of proportionally controlled air- 
craft. Aw = 0.25. Upper curve shows 6; lower shows #6. One 
div. = 0.5 sec 





Transient response of proportionally controlled air- 


Fic. 7. 
craft. \w = 0.5. Upper curve shows /; lower shows @. One 
div. = 0.5 sec 





Transient response of proportionally controlled air- 


Fic. 8. 
craft. Aw = 1.0. Upper curve shows 6; lower shows 6. One 
div. = 0.5 sec. 
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Part II—SupPeRSONIC AIRCRAFT 
(6) Equations of Motion 


Eqs. (6), (7), and (8) again describe the motion of the 
aircraft. However, an important simplification now 
can be made. It can be shown that the coefficient Z,, is 
approximately inversely proportional to the forward 
velocity or [Ep. Note: There is No Eq. (53)] 


Z.~ —g/U (54) 


from which it follows that for large forward velocities 
Z, becomes small. When evaluating this coefficient, 
an actual numerical example for a supersonic aircraft 
(which will be discussed later), the results indicate 
that Z, can be completely neglected. Furthermore, 
the assumptions of Eq. (9) can still be considered cor- 
rect. The important simplification arises now that the 
forward velocity perturbation u alone enters into .Eq 
(6), which means that one can work with Eqs. (7) and 
(8); and Eq. (6) serves only to determine u when de- 
sired. For convenience, rewrite Eqs. (7) and (8) as 


[p — (Z,/m)]a — [p — (g sin %/U)]@ = 0 (55) 


(o _ vp) 6 + - (—M,)a = a (56) 


and introduce the new aerodynamic coefficient 
, - 
2 = —Z,,/m (57) 


where 2 has the dimension of a frequency. The equa- 
tion of motion becomes now 


(p+ Qa — [p — (gsin&/U)]@=0 (58) 
(p? + 2Weuwwp)O + wy?a + Awy? = O (59) 


where Eq. (59) is identical to Eq. (22). Eq. (58) takes 
the place of Eq. (21), and an interpretation for this 
equation is in order. First, assume that the angle of 
attack of the plane is constant and soa = 0. (In the 
case of the conventional aircraft, this leads to the phu- 
goid theory.) Inthe present case 


[p — (g sin 0)/U)]@ = 0 (60) 
and 
0 = oye 8 0/0) (61a) 


In order to obtain a better picture, assume, for 
purposes of analysis, that at steady horizontal flight the 
attitude of the craft, , equals 3° and that suddenly the 
aircraft is made to assume an upward motion of 1° in- 
clination so that the attitude of the plane is made to 
change to 6; = 4°. (The angle of attack remains un- 
changed—that is, it stays zero.) Assuming U = 2,000 
ft. per sec., Eq. (61) becomes 


6 = b,e"" (61b) 
where 4) = 1,330 sec. = 22.2 min. From this, it 


follows that, for time intervals considered in these 


SUB, 19680 


dynamic stability studies, stays constant. In othe 
words, when the plane takes a direction different fro 
horizontal, it will keep this direction, unless, of course, 
new disturbance changes this direction. Thus, it cay 
be seen that the situation is different from the sy}. 
scenic case, where the phugoid oscillations develop, 
Actually, the phugoid oscillations also exist in th 
supersonic case, but the periods become so long that, 
for any practical purpose, the period can be considere; 
infinite. (If Z,, is not neglected, these long periods cap 
actually be computed.) 

It will be useful at this point to give a physical jp. 
terpretation for 2. The aircraft is assumed to be flying 
around a vertical loop (Fig. 9), and this is accomp. 
lished by keeping the angle of attack constant, and 
therefore a = zero. (This is equivalent to assuming 
that the lift force is much larger than the weight of the 
plane.) Neglecting g sin )/U, Eq. (58) becomes 


pd = Qa (62 


which means that Qa gives the radian-per-second tum 
of the plane. Thus, © 7s the angular rate of turn per 
unit of angle of attack change. For instance, if 2 =| 
rad. per sec. per rad., 1° of increase of angle of attack 
will make the aircraft turn 1° per sec. or 1 revolution 
in, 360 sec. For many purposes it will be useful to 
introduce the time constant given by the reciprocal of 
Q—i.e., 


Tt3 = 1/Q (63 


From Eqs. (58) and (59), the angle of attack a can be 
eliminated and the following equations obtained: 


E + 25 uWupP + Wy” = (g om A il ] 6 + 
p+oa 
New? 5 = 0 (64 


The last equation describes the exterior dynamics of the 
aircraft and takes the place of Eq. (23). 


Important differences between the low- and high- 
speed aircraft arise from the differences in the numerical 
values of the aerodynamic coefficients. As an example, 
the following coefficients might occur in the case of 4 
small supersonic high-altitude aircraft. 


Lift / 
Ys ‘ 
Lat = Const. 


7 — Weight 
“Suen on” 





Fic. 9. Supersonic aircraft flying along vertical loop. 
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Ww» = 20 rad. persec., % = 0 
t, = 0.06, 2 = 1.0 rad. per sec. per rad., A = 0.2 


When this is compared with the characteristics of the 
subsonic fighter, one is struck with the fact that wind- 
tunnel frequency w, is increased by a factor of 6, while 
the wind-tunnel damping {, is decreased to about one- 
tenth its subsonic value. This indicates that the super- 
sonic aircraft represents a fast oscillating system with 
little damping. To stabilize and steer such an aircraft 
will require a fast-acting well-matched servo system. 





(7) Uncontrolled Flight Fic. 10. Transient response of uncontrolled supersonic air- 
In this case, Eq (64) simplifies to craft. Upper curve shows 6; lower shows 6. One div. = 0.1 sec. 
,b — (g sin 4/U) ee, ¢ = 0.06 + 0.025 = 0.085 
E + 2 uWwP + ww" p a aa 6=0 (65) i : 7 : 
p+2 It can be seen that the frequency of oscillation in 


free flight is the same as the wind-tunnel frequency, 


wd while the damping ratio has a somewhat higher value, 

ore F 2 5» 2 though still small. As Eq. (66) is an equation for 

pe + 2 (‘- + “2 WP + Ww" 1+ 2% a pi =0 dé/dt, the rate of change of attitude of the aircraft, it is 
LW wo 





~ clear that the aircraft has no inherent tendency to re- 
turn to horizontal flight when disturbed. (The reason 
where it was assumed that g sin %)/U is negligible. This for this lies in the fact that the phugoid frequency is 
is a second-order equation for the rate of change of zero.) Assume, for instance, that a sudden angular 
attitude of the aircraft. The undamped frequency is velocity 6) is imparted to the aircraft. Fig. 10 shows 
the resulting oscillations of 6 and @. It can be seen 
/ ———___—— _ that the aircraft will achieve a steady-state sloping 

@ = Wy VI + 2fu(B/ww) (67) flight after the disturbance dies out. 


(66) 


given by 


and the damping ratio is given by 
(8) Proportional Control 


c - [Sw + (Q 2ww)] (Wy /w) (68) a . ~~ R ‘ 
To examine the possibility of using proportional con- 
For the supersonic aircraft considered in Section (6), trol, it will be assumed that 
w = wy = 20 rad. per sec. 5 = wd (69) 


Eq. (64) becomes, now, 


«i eieest 
G + 2uatap + wa 2 — Sem Osh) 5 dwn @=0 (70) 
p Q 
or 
? : Q F ft Q g sin % . al 
E +2 (‘ i 3 ) Wp” + Wy" (2 + 2_— + Mm) p + (rus saci U ) wef = 0 (71) 
aWry Wry 


For all practical cases g sin 0)/U < uw and so g sin 6)/U can be neglected. Eq. (70) is a third-order equation, 
and it can be studied with the aid of the cubic chart of the Appendix. The dimensionless parameters of Eqs. (131) 


and (132) are given by 


a = | > Au + 2fu (2/ww) =] * 1+ Aw (“*) (72) 
(Ap)? Q / (Ap)? Q 
of e 0 /I%..» »\ 78 9¢ \ */? 
Qs = 2[ Sw oth. 2ww) | (“*) ‘ais 4 = (“) (73) 
(Ap) “* Q (Ap) “* \Q 


where 
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: 
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Fic. 11. Transient response of proportionally controlled air- 





craft. \ = 2.0; w = 1.0. Lower curve shows vane motion. One 
div. = 0.1 sec 

Fic. 12. Transient response of proportionally controlled air- 
craft. A = 2.0; w = 2.0. Lower curve shows vane motion. 
One div. = 0.1 sec. 


is the approximate damping ratio for free flight. [Com- 
pare with Eq. (68), where w ~ w,.] First, a check of 
the stability of the system is made. According to 


Rouths'’ criteria, 


However, 


Wy “Wy b 
ajaa = 
Au 
(76) 
Q wy 
2ru 
20, 2 
= Se ME Re 
Au 


which shows that the system is stable indeed. 

The degree of stability of the system can be studied 
by examining the values of a; and a2. For most prac- 
tical cases, one is on the lower right-hand side of the 
cubic chart where w, is large. Using the appropriate 
approximate equations of the appendix, Eqs. (143) 
(145 = 
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A» f c 


2VAL Vit dwt (252/on) V1+ty 
(78) 


ae A, ] + Au + (252 Wy l 
‘an = = ~ l + ) 


Ao Aw Au 
(79) 
where 7; is the reciprocal of 2 [Eq. (63)]. 
Using the numerical values of the aircraft of Section 
(6), 


w, = 20V1+02n (80 

c. = 0.08/V 1 + 0.22 (81) 

T,, = [1 + (1/0.2p) ]73 (89) 

As an example, using uw = 1 yields w,, = 22 rad. per 
sec., ¢. = 0.073, T,, = 6 sec. This system is unsatis. 


factory; there is hardly any damping, and the elevators 
have to follow fast oscillations (3.5 cycles per sec 
But the time lag, 7,,,, is 6 sec., which indicates that, 
while the aircraft is oscillating with a high frequency, 
a great amount of time is required to steer the aircraft 
back to its course. Using a higher value of yu will in- 
crease 7°,,, but it will also decrease ¢,, making the sys- 
tem even more oscillatory. Figs. 11 and 12 show the 
transient response of the aircraft using » = 1.0 and 
2.0. Thus, it can be seen that proportional control is 
unsatisfactory. The following sections of the paper 
will present a more satisfactory type of control. 


(9) Fast Oscillations Versus Steering 
In order to understand better the control problem, 
some special considerations will be presented in this 
section. 
If the aircraft is suspended in a wind tunnel at its 
center of gravity, then 
a=06 (83) 
and the equation of motion [Eq. (59)] becomes 
(p? + 2Fuw@wP + @w")O + Aow"d = 0 (34 
Assuming that the elevators are controlled by a pro- 
portional servo—that is, 
6 = po (85) 
the equation of motion becomes now 
[p? + 2 uwwp + (1 + Au)on?]6 = 0 (00) 


The natural frequency of the system is given by 


w = wyeV1 + ru (87) 


When compared with Eq. (77), this shows that the fre- 
quency of fast oscillations equals the wind-tunnel fre- 
quency. On the other hand, the damping ratio of Eq. 


(S6) is 


& = Su/V 1+ Xp (90 
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which is smaller than the damping ratio of the free 
flight given by Eq. (78). 

Let us now turn our attention to steering. 
ing that a given elevator deflection instantaneously 
develops a certain angle of attack for the aircraft, the 


Assum- 


equation 
a = —ddb (89) 


holds. [As explained in connection with Eq. (27), this 
equation actually holds only for the static state.] 
From Eqs. (85) and (89), it follows that 


a = —pn0 : (90) 


and, from Eq. (58), 
[1 + (1 uA) |73p + 1} = () (91) 


where, again, g sin 0)/U is neglected. Eq. (91) is a 


first-order equation with the time constant 


Compared with Eq. (79), it can be seen that this time 
constant equals the time constant of the proportional 
servo-controlled aircraft. Summarizing, the dynamic 
parameters of the aircraft can be separated into two 
groups: 

(a) The fast oscillations correspond to the wind- 
tunnel condition (though the actual damping ratio is 
larger than the wind-tunnel damping ratio). 

(b) The steering motion corresponds to a highly 
simplified (nonoscillatory motion) characterized by 
Eq. (89), or the “immediate” angle of attack response 
case. 

The servo to be designed should perform two differ- 
(a) damp the fast oscillations and (b) per- 
Proportional con- 


ent tasks: 
form fast and effective steering. 
trol is unsatisfactory, since it does not damp the fast 
oscillations and does not steer fast enough. Further- 


more, Eqs. (78) and (79) show that speeding up steering 


or 


j l + 26 wl wre 


p +. 





p? 


results. 

In order to study the behavior of this cubic equation 
with the aid of the cubic chart, a; and a» of Eqs. (131) 
and (132) are computed: 


260 -+- ( T; + Aull, wr - 
a ae W374 i (97) 
(1 + Ap) (Tw) ”’ 
9+ (T . 
a = l + 26 w( Tw) (98) 


(1 + Ap) “(Tw)”? 
It is important to notice that 47, enters only into the 
formula of a. This means that, when varying u7); 
one moves on the cubic chart on a horizontal line, and 
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speeds up the steering. 
will be arrived at by separating the two problems 
that is, first by developing a servo system that damps 


= + (2¢. + (T, = 
T; Tl; 
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(increasing yu) will result in decreasing the damping of 


the fast oscillations. 


in the remainder of this paper, a control system will 


be developed which damps the fast oscillations and 


The proper interior dynamics 


the oscillations in the wind tunnel [Eq. (86)]; then by 


developing a servo system that steers speedily [Eq. 


(91)]; and, finally, by showing that a combination of the 


systems will result in a satisfactory performance for the 


complete aircraft. 


(10) Theory of Fast Oscillations 


A study will be made of the appropriate control func- 
tion to stabilize the motion described by Eq. (S4). 
Because of the high frequencies involved, the time 
lag of the servo will be important. In order to simplify 
matters, it will be assumed that the relationship be- 


tween signal and vane motion is given by 
(T.p + 1)6 = wo (93) 


This means that the servo is assumed to have a simple 
exponential time lag, 7). Furthermore, it will be as- 
sumed that o, the signal to the autopilot, is a combina- 
tion of the attitude and the derivative of the attitude 


of the missile; thus, 
6 = [a + uTap)/(L + Tip) | (94) 


We will try to determine » and u7,* for optimum per- 
formance for any given fixed time lag, 7). Substitut- 


ing Eq. (94) into Eq. (S4), 


(»* + etoop + wet + Mune? 2 rae) @=0 (95) 
1+ Tp 


* For the present, «7°, should be considered as a single param- 
eter with the dimension of time. Later, it will be factored into a 


dimensionless coefficient and a time factor 


AuT, wre lou p+ (1+ dp) ot =u f (96) 


1 


therefore, for any value of a2, a value of wl, can be 
determined to give the maximum possible damping 
ratio. This maximum damping ratio is given by Eq. 


(147) of the Appendix; thus, 


or T.. 1 
= (0.192 LL + 2f007 we] (99) 


T wu V1 + dB 


fopt. = 0.192 az 


The last equation shows that @ = 0 gives more damping 
than any positive @ Physically speaking, this means 
that the servo signal should not contain any propor- 


tional term but only a derivative term. Eq. (99) be- 


comes now 





~~ 
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5, = 0.198 (1+26,T, wn)” 
. =0. 
Wn 


L 





6 7 86 9 


Ty Wp 





Fic. 13. Optimum damping ratio vs. inherent damping ratio and Ty. 


Sop. = 0.192 -- (100) 


T IW 


or 7 a/s 
(1 + 2ST ww) ' 


Fig. 13 shows a plot of the last equation. For small 
values of ¢,, small values of Tw, are required in order 
to leave fair control damping ratios. For the aircraft 
considered before, ¢, = 0.08. In order to have, say, 
¢. = 0.25, one must have Tw, = 1.0. 
T, = 0.05 sec., which shows that an extremely fast servo 
is required to damp the aircraft satisfactorily. 

It might be of some interest to compare these co- 


For aw, = 20, 


efficients with the conventional aircraft treated in 
Part I of the paper. One has ¢, = 0.685 and w, = 
3.03. Fig. 13 shows that, whatever value of the time 
lag of 7; is assumed, the control damping is above 
0.70. Using, for instance, Tw, = 2 or 7; = 0.663 sec. 
will result in ¢, = 0.71. Stabilization of an ordinary 
aircraft in the wind tunnel does not require a fast servo, 
and the speed of the servo is determined by the desired 
response time to return to equilibrium. 


(11) Theory of Steering 


Eqs. (79) and (91) show that, in order to speed up 
steering, a large value of u must be used, where yp is the 
ratio of elevator deflection to attitude in static case. 


However, it is known from Eqs. (78) and (91) that in- 
creasing u decreases the damping of the fast oscillations. 
In order to avoid this bad interaction between the two 
types of motion, a proportional control is introduced— 
that is, 

6 = [u/(1 + Tp)]o (101 
where 7 is a time lag. The value of this lag should be 
large enough so the elevator does not follow the fast 
oscillations of the aircraft but not so large that efficient 
steering cannot still be effected. As a try for the value 
of T, use 


T = 1/2 = 73 (102 


The equation of motion, Eq. (61), becomes now 
(»* + ennap + AE _ 4 se) 6=0 (103 
Q + p Q+ p 
or 
[p? + (2Fnww Q)p? + (ww? + Welw) Pp + Aporw?Q) X 
¢=0 (104 


The performance of this system can be studied with the 


aid of the cubic chart. Thus, 
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Q wy \'”? Wy \ 
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- ( ' Wy) \ApQ AuQ ) 
9) 7 Vs 
x= (04 3)" 
WJ \ApQ 


since, for most cases, 2fm(Q/w») << 1. 

For most practical cases, one is on the lower right- 
hand side of the cubic chart, where w, is large, and Eqs. 
(143)-(145) of the Appendix can be used. Thus, with 
the aid of Eq. (74) and (104) and assuming that 2¢,,2 + 


Wy < i 


(106) 


+ = As ‘2A, => eo + (2/2) = a (107) 


and 
Wn, = Ay = Wy V1 + (2¢,,0 Wy) =< Wy (108) 
1 2¢,,2/ a.» 
lw = A, Ao = [PE ehaie 1 = 73 (109) 
Au Q Au 


By comparing these equations with Eqs. (74), (75), and 
(76), it is seen that the delayed proportional control is 
indeed superior to the straight proportional control, 
since the former does not decrease the damping ratio and 
results in a small-time exponent for steering, depend- 


G + 25 Ww + ( 


This gives the following fourth-order equation: 


{ 


l 


9) _ L@) 

“Ww wrt 

a aes a 
T; 


In order to study the performance of the aircraft, the 
roots of Eq. (112) will have to be determined. The 
procedure to be followed for the specification of the 
servo system can be illustrated by considering the air- 
craft discussed in Section (6) of this paper. 

It is arbitrarily specified that the damping ratio of 
the system (that is, the damping ratio of the oscillation 
component) should be of the order of 0.20. More spe- 
cifically, it will be required that the “wind-tunnel” 
damping of Eq. (100) should be 0.20. Then ¢, = 0.06, 
and therefore, from Fig. 13, Tw, = 1.3. As wy, = 
20.0 rad. per sec., the time lag of the servo is 7, = 
0.065 see. This figure shows that a fast servo is re- 
quired to properly stabilize the system. Now, deter- 
mine the value of 47, of Eq. (110). With the aid of 
Eq. (98) (@ = 0), a2 = 0.975. According to the cubic 
chart for this value of as, a; = 3.0 in order to satisfy the 
assumed ¢ = 0.20. Solve Eq. (97) for the unknown 
xT, and obtain u7, = 0.5 sec. 
of Eq. (110) is given by 


Au + tsp 
1 + 73p 
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ing on the value of Au. [Note of caution: This last re- 
sult depends on Eqs. (107)-(109), which hold when w, 
of the cubic chart is large enough. When w, is not too 
large, this is usually the case. However, the upper 
limit of Au, and therefore the minimum of 7,,,, is given 
only by the requirement that w, must be large enough. } 

To summarize, delayed proportional control gives 
good steering and does not affect the fast oscillations. 
The derivative control described in Section (10) im- 
proves the damping of the fast oscillations and does not 
affect steering, since there is no term proportional to 
the attitude of the missile involved. A combination of 
these two types of controls should result in a satisfactory 
control system for the supersonic aircraft. 


(12) Tuned Computer-Network-Controlled System 


In accordance with the results obtained, the follow- 
ing control system is proposed for the aircraft: 


6 = [uTap/(1 + Tip)]6 + [u/(1 + rsp)]@ (110) 


where the first term is the ‘‘damping’’ term and the 
second is the “‘steering’’ term. When substituting 
Eq. (110) into the equation of motion, Eq. (64), the 
following equation is obtained: 


J ps rm (2500 - Q 4 - )o* + oo Au ntti + 26 Wy (2 + = )| p? + 
; U 


Mita ot 0 = 0 (111) 

1+ Tip 

Wy? a1+TM% , | Wy? Q) ;' 

— + Ap —-_ o, 72 +X ——>§@=0 (112 

r, ld T, w p vad T, { ( ) 
5 = [0.5p/(1 + 0.065p) ]@ (113) 


It will be assumed for the time being that steering 
is completely independent of the fast oscillations. 
Thus, as far as steering is concerned, it will be assumed 
that the wind-tunnel damping ratio is not 0.06 but that 
fo» = 0.20, the damping coefficient accomplished by 
using a control given by Eq. (113). Using this value 
of ¢ in Eqs. (105) and (106), then, 


a, = 1.02 (20/Au)*”’ (114) 


az = 0.45 (20/du)'” (115) 


A selection of wu is made such that the system has a 
(In effect, 
Selecting 


damping ratio near to the original 0.20. 
select » such that w, remains to be “‘large.’’) 
uw = 10, ay = 4.75, ag = 0.97, and, from the cubic 


chart, ¢ = 0.18. Thus, the proposed steering system 


Thus, the first term * ” 


6 = [10/(1 + p)]@ (116) 
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According to Eq. (109), the time constant of the steer- 
ing system will be approximately 7,,, = 0.5 sec. 
The complete control equation is given now by Eq. 
(110) with the following coefficients: 
T, 0.065 sec., yw = 10) 
73 = 1.0 sec. (117) 
T 0.05 sec. 


Using these values, the coefficients of Eq. (112) can be 


p' + 18.863 + 1,019p? + 7,600p + 12,300 = [p + (1 


From this it can be seen that the damping ratio of 
second-order component is indeed 0.18 and that the 
undamped natural frequency is 30.6 rad. per sec. The 
fourth-order equation has two other real roots, giving 
the exponential time constants 0.435 and 0.175 sec. 
(The first of these corresponds to the value of Pays esti- 
mated to be 0.5 sec.) 

Fig. 14 shows the transient response of the aircraft. 
It can be seen that good damping and fast steering are 
achieved but that a fast-acting servo is required. 


(13) Computer Design 


A few words may be said about the design of a system 
realizing the control equation [Eq. (110)]. Either a me- 





Fic. 14. Transient response of tuned computer controlled 
supersonic aircraft. 











Fic. 15. Computer circuit 
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determined. In order to check the validity of the as. 
sumption that the control systems are really inde. 
pendent of each other, the roots of Eq. (113) are to be 
determined. A numerical calculation gives* 


* Only a single real root of the fourth-order equation has to be 
determined, and, as a first approximation, 1 T,, can be used 
When this root is determined and factored out, the roots of the 
resulting cubic equation can be determined with the aid of the 
cubic chart. 


0.435)] [p + (1/0.175)] X 
(p? + 2 X 0.18 X 30.66 + 30.62) (118 


chanical or ehevtoonie computer could be used. Assume 
that the sensitive element of the autopilot gives a vol- 
tage proportional to the attitude of the aircraft—that 
is, 


V = pw (119) 
and assume the design to be an electronic computer. 


When comparing Eq. (120) with Eqs. (93) and (110), 
one obtains 


(120 


c= (ap + 


1 + aa , 
Mip2 


1+ T3P 
which is the operation to be performed by the autopilot 
computer. 

Fig. 15 shows the schematic design of such a com- 
puter. An easy circuit analysis gives 


Ey = }[1+ RC,)/{1 + (Ri + RC IIE, (121 

By selecting 
RC =T, (122 
(Ri + RC 


(128 


II 


Ts 
3 


the second term of Eq. (120) is obtained. The first 
term of the equation (the derivative term) can be real- 
ized only approximately. 
Taking 
the voltage across the resistor is 
Ey = RC[p/(1 + RC,)JE, (125 
When 
RG, <1 (126 
the circuit gives a derivative term. By using two 
different circuits of the type of Fig. 15, amplifying the 
signals, and adding the voltages, the operation of Eq. 
(121) can be performed. 
The reader who is interested in the actual design of 
the computer is referred to the appropriate literature on 


computers. * 
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CONCLUSIONS 


It will be of some interest to present a few general 
words on the method of preparation of an analysis and 
synthesis for an autopilot-controlled aircraft. 

“When the engineer is confronted with the problem, 
the first thing is to prepare a block diagram. How com- 
plete this block diagram should be is difficult to specify 
and somewhat up to the individual taste. Generally, 
the more complete the diagram the better. 

The next step is to obtain a thorough understanding 
of each box of the diagram. It should be clear what 
each box means and what function is performed by each 
box. Complete understanding of each box is obtained 
only after the transfer function (or the system of dif- 
ferential equations connecting input and output) is 
obtained. Combination of these differential equations 
forms the mathematics of the system. 

Usually, the system of differential equations obtained 
by this fashion is complicated and of high order. It 
would be difficult to start the analysis on the basis of 
such a complicated system of differential equations. 
For this reason, it is of extreme importance to try to 
simplify the system and thereby analyze by parts. It 
is difficult to give specific instructions on this matter of 
simplification, but, generally speaking, one is inclined 
to tend to neglect too few things at the beginning. 
Some simplification might look too absurd. However, 
after a simplified analysis is made, it might turn out 
that much valuable information is obtained at little 
expense in time and effort. Furthermore, information 
can be obtained on the effect of a single box of the block 
diagram, and often the mathematics of a single box can 
be replaced by simpler mathematics, making the analy- 
sis of the whole system much easier. 

As an example, we can point out that, in the analysis 
presented in this paper, it was assumed that the auto- 
pilot itself was governed by a simple time lag. A 
much more complicated analysis (not presented in the 
paper) showed that the most important features of the 
design could be evaluated on the basis of this simpli- 
fied theory. 

After a simplified analysis is completed, one endeav- 
ors to refine the analysis by including more and more 
of the actual physical properties of the system. How- 
ever, one soon arrives at a dead end. The equations 
become of extremely high order and nonlinear. 

Nonlinearities in autopilots arise from many sources. 
For instance, the elevators are always restricted in their 
motion by stops. Hydraulic systems have finite size 
ports and other nonlinearities, electric motors have 
speed limitations, and vacuum tube outputs are limited. 
These nonlinearities have a strong effect on the per- 
formance of the system and put a limitation on the 
value of the linear analysis. In particular, a refined 
linear analysis might be illusory, since it might happen 
that the effects due to the refinement of the analysis 
are completely overshadowed by the nonlinear effect. 
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There does not seem to be available at the present 
time a general theory of nonlinear systems, permitting 
the analysis of nonlinear autopilots. On the other 
hand, computing machines and, in particular, analog 
computers are available or can be built which can 
handle sufficiently general nonlinear autopilots. There- 
fore, in most cases, the final analysis of an autopilot 
can be performed with most economy on a computing 
machine. 

It appears that the mathematical analysis serves to 
obtain an insight and a quantitative understanding of 
the problem. When a complete analysis is required, a 
computing machine analysis should complete the paper 
work, 

In many cases, the designer of the autopilot will not 
require a complete analysis. After the simplified paper 
analysis gives him an insight into the problem, he will 
proceed to build the autopilot and the computer with 
adjustable elements. Experience shows that, however 
refined the analysis is, the final adjustment of param- 
eters and the checking out of the system must be done 
on so-called flight simulators or in actual flight testing. 

As a summary, the following work schedule can be 
set up: 

(1) Block diagram. 

(2) Analysis of blocks of block diagram (component 
analysis). 

(3) Mathematical analysis of simplified system. 


(4) Computing machine analysis. 
APPENDIX 


(A) The Cubic Equation 
This equation is of the form 
(p* + Aop* + Aip + Ao)é@ = 0 (127) 


There is always one real root that can be factored 


out. Thus, 


p*? + Asp? + Aip + Ao = [p + (1/T,,)] X 


(p? + 2f3wn,P + wn) (128) 


where 7°,, is the natural time lag of the system, w,, is 
the (undamped) natural frequency of the second-order 
component, and ¢3 is the dimensionless damping ratio 
of the second-order component. 

The performance of the cubic equation was studied 
by Y. F. Lin with the aid of the cubic chart (Chart 1).° 
By introducing the dimensionless time 


= VAg (129) 
Eq. (124) transforms into 
ae 
where 
a, = A;/Ao”’ (131) 








> 


Ce ae 
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and 


a = A2/Ay”* (132) 


The characteristic equation corresponding to Eq. (130) 
is 


+ adr + art+1=0 (133) 


and it can be factored into a first- and second-order 
factor* 


3+ awd? + ard + 1=[A+4+ (1/w,?)] X 


(A2 + 2¢30,A + w,2) (134) 


The cubic chart gives w, and {3 as functions of a, and 
a. This chart was prepared by using the formulas de- 
rived from Eq. (134): 


a = w,? + (263/e,) (135) 
a, = (1/w,") + 2630, (136) 


A further important relationship can be obtained by 
comparing Eqs. (127) and (134): 


Wn, = at Me (137) 
Tu, = 0,3/V Ap (138) 
Assuming 
w, > 1 (139) 
the above equations simplify considerably. In_par- 
ticular, Eqs. (135) and (136) become 
a, ~ wo,” (140) 
an ~ 263w, (141) 
bs ~ an/2V on (142) 


Furthermore, with the aid of Eqs. (131) and (132), it 
follows that 


wn, ~ VA (143) 
Tn, ~ A:i/Ao (144) 
bs ~ A,/2V A, (145) 


Thus, it can be seen that, for the w, ‘‘large’’ region, the 
characteristics of the cubic equation can be deter- 
mined without the aid of the cubic chart. When these 
approximate formulas are used, an a posteriori verifi- 
cation is necessary to show that w, is indeed sufficiently 
large. 


Another important formula for the cubic equation is 
the maximum possible damping ratio for a given az. 
This is given by the minima of the constant ¢ curves 
or 


* 


wr is a dimensionless quantity. The notation is taken from 


Lin’s paper.® 
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(dae/doai)y = cons, = 0 (146) 


A simple calculation leads to the following results: 


S masz. = (a2 3) Pa 0.192a2 ¥ (147) 


and 


a, = 3/ae + 2/9ae' 


(148) 
(B) Statically Unstable Aircraft 


When the center of the aerodynamic forces is ahead 
of the center of gravity, the aircraft is statically yp. 
stable. If the control surfaces are locked, the air- 
craft will tumble. It is an interesting problem to ap. 
alyze the stability of the aircraft when controlled by an 
autopilot. , 

In order to understand the problem, we will assume 
that the aircraft is suspended in a wind tunnel at the 
center of gravity. The dynamics of the system can 
be described by 


(p? + 2owp — w?)? + ws = 0 (149 


This equation replaces Eq. (26), which holds for the 
statically stable case. 

When the control surfaces are locked, Eq. (149) be- 
comes 


(p? + 2p — w?)d = 0 (150) 


which shows that the motion is of the exponential type 
with a positive time coefficient as the exponent. A 
physical interpretation can be given to Eq. (150) by 
considering, as a model, the behavior of an inverted 
pendulum suspended below the center of gravity. 

The simplest control system to be investigated is as- 
sumed to be proportional or 


5 = 70 (151 
Eq. (149) becomes 
[p? + 2p + (Ad — 1)w7]@ = 0 (152 
The condition for stability is 
z>1/r (153) 


which simply means that the autopilot ‘‘spring-constant” 
must be larger than the destabilizing aerodynamic “‘spring- 
constant.” 

Let us now investigate the effect of a time lag. For 
this purpose a control equation similar to Eq. (94) is 
assumed : 


6 = [(@ + wTsp)/(1 + Tip)] 6 (154) 
Eq. (149) becomes 
Agw* + netee't) ¢@=0 


(»* + 2wp — w* + — 
1+ 7p 


(155) 


or 
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Jp 1+ 2foT'w p+: 
\! 3 5 


This third-order equation can be studied with the aid 


of the cubic charts. Using Eqs. (131) and (132), 


2¢ ul, — T)) a 

‘aise fo + ( caw SS. w (157a) 

(Ag — 1)7*°(Tw)” 

. oe.T 

ac= : + 2foT wo = (157b) 

(Ag— 1 ) “ # 10) die 

The best possible damping ratio for Eq. (144), 

¢ = 0.192 (1 + 2607 w)'*/TwV/dg — 1 (158) 


S opt 
Fig. 13 shows a plot of this function without the factor 
V \a — 1 in the denonnnator. 

At a first glance, it would appear from Eq. (158) that 
good damping can be obtained by making Aj slightly 
larger than unity. However, a more detailed analysis 
shows that minute changes in the parameters can com- 
pletely destabilize the system for this case. As for 
any actual physical system the parameters will vary 
between certain limits, one concludes that Ag must be 
sufficiently larger than unity. For purposes of argu- 
ment, it can be assumed that A\g = 2. 


Eq. (158) becomes 


(159) 


(1 + 26oT w) 


Cont, = 0.192 
Tw 


which is plotted on Fig. 13. Two cases are to be dis- 
tinguished. 

(a) Fair value of fo, say, {9 > 0.5. Conventional 
subsonic aircraft falls into this group. 
that there is no special requirement on the value of 
Tw. A statically unstable conventional aircraft can be 
controlled by an autopilot without demanding unusually 


high speed of autopilot action. 


> 


Fig. 13 shows 


(b) Small value of ¢>. Guided missiles under cer- 
tain conditions fall into this group. 
in this case 7 must be small. When considering ac- 
tual numerical cases, the following conclusion can be 
For a supersonic aircraft, a small amount of 
Larger static in- 


» 


Fig. 13 shows that 


drawn: 
static instability can be compensated. 
stability can be compensated only by a fast-acting auto- 
pilot. 


(C) Aircraft Controlled by Rotating Wing 


Some of the peculiarities of aircraft control arise 
from the fact that the elevator deflection does not di- 
rectly change the lift forces of the aircraft. This situa- 
tion can be remedied by mounting the wings on a shaft 
and changing the incidence of the wing itself. It will 
be of interest to establish some of the basic aspects, 


of this type of control. Only the supersonic case will 


-CONTROLLED 


[2f0 + nu Ts N xs 
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Polo p+ (a (156) 


be investigated. For sake of simplicity, it will be 


assumed that the tail configuration is fixed. 


Denoting by 6 the change in angle of attack of the 
wing, Eqs. (6), (7), and (8) remain valid. However, the 
normal force coefficient Z; is now important. In order 
to simplify matters, it will be assumed, when the wing 
is rotated, that the body axis remains unchanged or that 
no moment is exerted by the wing. This shows up in 


Eq. (8) by neglecting 1/;. Eqs. (58) and (59) now be- 


come 
(p+ Qa — (» _§ sins) 6+2;6=0 (160) 
(pe + 2FywoP)O + wy?a = 0 (161) 
where 
Q; = —Z;/mU (162) 


The ‘‘wind-tunnel’’ equations follow from Eq. (161). 


(p? + Wuwwp + wy?)0 = 0 (163) 

Because the wing exerts no moment, the wind-tunnel 
characteristics are independent of the autopilot. In 
particular, if there is not sufficient aerodynamic damp- 
ing available, the autopilot cannot remedy the situa- 


tion. 


From Eqs. (160) and (161), @ can be eliminated. 


Thus, 
ju aes Pp — (g sin @/U) 
4 2 uw dn 2 wT = 
| P+ 2uneb + « +o x 
Q; 
6+ w,? 6=0 (164) 
i Q 


This last equation describes the exterior dynamics of the 
system. When compared with Eq. (64), Eq. (164) can 
be interpreted as a The 
“theory of steering’ of Section (11) applies directly to 


“delayed” type of control. 


this case by putting 


(165) 


Using this interpretation of our “theory of steering,” it 
can be shown that movable wing will give good steering 
but no increase in damping ratio. When increasing the 
wing effectiveness Q;, the steering becomes faster, but 
the pitch damping coeflicient decreases. Whether wing 
control can be made satisfactory will mainly depend on 
the availability of good aerodynamic pitch damping. 


(D) Analysis of Aerodynamic Coefficients 


The essential difference between subsonic and super- 
sonic aircraft lies in the differences in the values of the 
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aerodynamic coefficients and, in particular, in the 
values of the wind-tunnel frequency w, and damping 
ratio ¢,. It will be of some use to explain this situa- 
tion, at least qualitatively. 

Suppose the aircraft is suspended at the center of 
gravity in the wind tunnel. The equation of motion is 
given by 


I (d?a/dt?) + D(da/dt) + Ma = 0 (166) 


where A/a is the aerodynamic moment and D(da/df) 
is the aerodynamic damping. M/ and D are given 
by 


Ma = (dC,/da)aAl(pU?/2) 
D(da/dt) = C(da/dt) Al?pU 


(167) 
(168) 


where dC,/da is the lift coefficient slope, A is the char- 
acteristic surface area, / is the distance between the 
center of gravity and center of pressure, / is a character- 
istic ‘distance, C is a dimensionless coefficient, p is the 
air density, and LU is the forward velocity. The wind- 
tunnel frequency and damping ratio is given by 


I Pa IC. l Al’ ; : 
On = . l Vp Pp L ( ) ~ { (169) 
V2 Via WIV wie 
— da / - 
Sn» = C V p b : ; ~ Vp (170) 
VO dC. VIAI 


Consider the variations in w, and ¢, for a given air- 
craft. Fora supersonic plane, there is a great variation 
in forward velocity, U, from the take-off to maximum 
velocity. dC,/da is a function of Mach Number, but 
the overall change is moderate. The distance between 
center of pressure and center of gravity varies con- 


siderably, since the center of pressure position is a 
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function of Mach Number and the center of gravity 


If the air. 
craft rises to high altitude, changes in air density be- 


position varies when the fuel is being used. 
come important. Also, the inertia / changes when the 
fuel is used up. Qualitatively, the following obserya. 
tions can be made: 

(1) High velocity results in high frequency. 

(2) Low air density results in low frequency. 

(3) Low air density results in low damping ratio, 

Finally, consider the effect of the overall size of the 
aircraft. 
tional to the fifth power of the dimension; 


Roughly speaking, the inertia / is propor. 
the area 
A is proportional to the square. 
Al*/I and /?/+/JAI are independent of size. 
Eqs. (169) and (170), it follows that: 

(4) Small aircraft results in high frequency 


Thus, approximately 
From 


(5) Damping ratio is independent of size. 

The least favorable condition arises when a small 
aircraft takes off from ground and rises to high alti- 
tude and high velocity, since in this case the frequency 
will rise to a high value and the damping ratio will be- 
come small. 
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Bifurcation Criterion and Plastic Buckling of 
Plates and Columns’ 


C. E. PEARSON?t 


Brown Unwwersity 


ABSTRACT 


If a bifurcation of equilibrium paths occurs as the loading is 
varied on a conservative system, the preferred path is usually 
obtained by use of a stability criterion. However, even if the 
definition of stability is modified so as no longer to depend on the 
existence of a potential energy function, a stability criterion may 
not be useful if the system is nonconservative. The choice of a 
different type of criterion is discussed for the particular case of 
plastic buckling of plates and columns, and detailed calculations 
are given for the simplified column model introduced by Shanley. 

To find the deflection of a real column with arbitrary cross 
section in the neighborhood of the tangent modulus load, numeri- 
cal methods must be used. Two such methods are outlined: 
One consists of a sequence of approximations to the simultaneous 
solution of two nonlinear integral equations, and the other con- 
sists of a sequence of graphical integrations based on the stress- 
strain curve of the material. Under certain conditions, a portion 
of the solution may be found in analytic form. 

The problem of plastic buckling of plates under compressive 
edge thrust has been most rigorously treated by Handelman and 
Prager. If their calculations are modified by the assumption of 
initial loading (in analogy with the tangent modulus theory for 
columns), then their results may be brought into fair agreement 
with recent experimental work by Pride and Heimerl. The as- 
sumption of initial loading may be verified a posteriori. 


STABILITY OF NONCONSERVATIVE SYSTEMS 


Fe THE ENGINEERING POINT OF VIEW, a structure 
is stable if the disturbances to which it is subject 
are not sufficient to produce excessive displacements. 
However, to permit mathematical treatment of stabil- 
ity, it is customary to idealize in two ways: The struc- 
ture itself is mechanically idealized (e.g., friction in 
pin joints, finiteness of some dimensions, etc., are often 
neglected); and only such disturbances are considered 
as are certain to be applied to any structure because of 
the vagaries of nature—.e., infinitesimal forces. 

If the idealized system is conservative, then the state- 
ment that infinitesimal forces will produce no finite 
deformation is equivalent to the statement that the 
potential energy of the position in question is a relative 
minimum. This latter statement is essentially the 
classical definition of stability; if the idealized system 
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* This paper is a condensed version of a thesis submitted in par- 
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is nonconservative, the definition of stability must be 
given with respect to infinitesimal disturbances instead 
of a potential energy function. The following definition 
is valid for both conservative and nonconservative 
systems: 

“A mechanically-idealized system will be said to be 
unstable if there exists a small distance ¢ such that for 
any 6 > 0 it is possible to find a force distribution of in- 
tensity less than 6 which will produce an eventual dis- 
placement of the system in which at least one material 
point is distant « from its original position. A stable 
system is one that is not unstable.”’ 

Consider now an arbitrary mechanical system, either 
conservative or nonconservative when idealized. If 
the loading is varied, there can be no ambiguity in the 
equilibrium path followed by the real system, but such 
ambiguities may occur in the analysis of the idealized 
system. In the conservative case, the preferred path 
may usually be found by considering the stabilities of 
the alternative paths, although, because the exact dif- 
ferences between the real and idealized systems are not 
known, it is not possible to prove rigorously that the 
real system will follow this preferred path. Of course, 
the less the idealization, the greater the likelihood that 
this will be the case. However, stability need not be 
the only criterion; it might even happen that several 
of the possible paths are stable, and in this eventuality 
an additional criterion must be sought. For a non- 
conservative system, the same remarks will apply, pro- 
vided, of course, that the modified definition of stability 
is used. 

In the problem of interest here 
there is no difficulty if the buckling 


viz., the buckling of 
columns and plates 
load occurs within the elastic range of stress, for a 
straightforward calculation shows that only the buckled 
state is stable. If the stress is beyond the elastic limit, 
not only is the stability criterion difficult to apply but 
it seems likely that both the buckled and the unde- 
formed states are stable; a different criterion is there- 
fore required. It is well known that it is often useful to 
examine the deformed state of a plate or column; if an 
initial deflection is introduced into the analysis, it will 
be found that the equilibrium path is always unique 
and approaches the buckling path as a limit when the 
initial deflection approaches zero. In other words, the 
convergence is nonuniform; the straight-line equilib- 
rium position is obtained only as a result of zero initial 


deflection and is not approached by any sequence of 
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i Fig. 3 Similar Triangles 
equilibrium paths corresponding to diminishing initial : 1 (EA E,7A 
deflection. This result is the basis for the use of initial aniltelias a\ 1 + - ] Al + 
deflection as a criterion. Of course, other criteria are 1 


possible—e.g., an eccentricity in the application of the 
load—but calculations are not so simple. 

To illustrate these remarks, the simple column model 
of Shanley' will be analyzed; it may be remarked that 
this column follows the behavior of more conventional 
columns more closely than might be expected. 


THE SHANLEY MOpDEL 


Consider the column shown in Fig. 1; it consists of 
two rigid sections, each of length L/2, joined by means 
of two plastic rods of lengths /; (inner) and /, (outer), 
It is assumed that at a load P, (less than the tangent 
modulus load P;) the deflection of the center is 6,, and 
that at a future load P the deflection is 6. Corre- 
sponding to this new deflection, the lengths of the plas- 
tic rods will be /; + Al;, 1, + Al,, respectively. The 
stress-strain curve for the material is that shown in Fig. 
2; the elastic modulus is /, and the tangent modulus 
(assumed constant) is £7. The width of the column is 
h and the area of each plastic bar is A. 

If the unstressed length of each plastic rod is /, then 
the increase in compressive force corresponding to a 
change of length Al is 


( EA. Bet) ud a 
2 I I , 


Now (d/dP) (Al;) is always < 0, and, if 6, is suffi- 
ciently small, then (d/dP) (Al,) will be < 0 for the ini- 
tial part of the motion. The equilibrium equations for 
forces and moments will read 


P— P, = — (E7A/l) (Al, + Al) (2 
and 
Pé — P,6, = —(E7A/l) (Al, — Al,) (h/2) (8 
From similar triangles (Fig. 3), 
Al, — Al, = (4h/L) (6 — 6,) (4 
Eqs. (3) and (4) may be solved for 6 to give 
6 = 6,(Pr — P,)/(Pr — P) (9) 
where P 7, the tangent modulus load is given by 
P, = 2E,Ah*/Ll (6) 


This solution must satisfy (d/dP) (2L cos a + Aly + 
Al,) < 0 (in order that the increment in load do post- 
tive work), (d/dP) (Al;) < 0, and (d/dP) (Al,) < 9 
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PLASTIC BUCKLING OF 
(see previous assumption). It is easily checked that 
the first two conditions are always fulfilled but that 
(j/dP) (Al,) is < Oonly if Pis < P, where 

P, = Pr — V(26,/h)Pr(Pr — P,) (7) 
This last equation is obtained by solving Eqs. (2), (4) 
5) for Al, in terms of P and by setting the deriva- 


and ( 
Note that 


tive of Al, with respect to P equal to zero. 
P.< Pr. 
For P > P., (d/dP) (Al,) will be > 0, and the equa- 


tions of force and moment equilibrium will now read 


P — P, = —(E7A/l) Al,’ — (EA/I) Al,’ (8) 


and 
Pi — P.6. = (Ah/2l) (—E,7Al,’ + E Al’) (9) 


where 6, is the value of 6 corresponding to P, and the 
primes mean that the changes in length are now to be 


measured from the load P,. Eq. (4) becomes 


Al,’ — Al,’ = (4h/L) (6 — 6,) (10) 


0 


Eqs. (8), (9), and (10) may be solved for 6 to give 


rs - (Fes a) +5 ( — “)(72=*F) (11) 
a °e Pr = P Z Pr = Pr 


where Pp, the reduced modulus load, is given by 
Pp = P7(2E/(E + E7)] (12) 
Using 6, = 6,(Pr — P.)/(Pr — P.) in Eq. (11) gives 


_. (Pr—P,) (Pr— P.) , h (? ™ _) 
Cn DI a PO SN 


=) 
(13 
( P, 
Note that P, < P.< Pr< Pp. 


The solution [Eq. (13) | must satisfy, as before, (d/dP) 
(Lcos a+ Al,’ + Al,’) < O and (d/dP) (Al,’) < 0; in 
addition, (d/dP) (Al,’) must be > 0. It may be 
checked that this is the case for P > P,. Because of the 
signs of Al, and Al,’, Eq. (5) is valid only for P < P.and 
Eq. (13) is valid only for P 2 P,; it is therefore clear 
that Eqs. (5)-(13) are unique. One special case re- 
mains—viz., 6, = 0. In this case the previous for- 
mulas show that for P < P, the solution is uniquely 
6 = Oand that at Pa bifurcation of equilibrium occurs 
with the two possible solutions 6 = 0 and 


sl h @ me *2)(7* i fe 
°~ 9\P,— P/\” Pr 


In fact, a bifurcation of equilibrium occurs at every 
position on the 6 = 0 path for which the load lies be- 
tween Pp and Pp; if the load corresponding to the bi- 
furcation point is P,, then the possible paths are 


6=Oand 
= (P= Pa) Pe = Pe) 
— 2\Pp — P P, 
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Now the previous definition of stability is difficult to 
apply; furthermore, it seems likely that both paths in 
the 6, = O case are stable (note that, for the 6 = 0 
path, no adjacent equilibrium position exists under 
constant axial load). However, irrespective of whether 
the paths are stable or not, the fact that only one of the 
equilibrium paths is approached as a limit of the 6, # 
0 solution is sufficient to indicate strongly that buckling 
will occur in the neighborhood of the tangent modulus 
load. This nonuniformity of convergence will be 
adopted as a bifurcation criterion rather than one 
based on an extension of the idea of stability (for an 
attempt along these latter lines see reference 2). Refer- 
ence is made to the discussion of the preceding section; 
see also the discussion of von Karman following refer- 
ence 1. 

In setting up the above equations, the moment arm 
of the forces in the plastic rods has been taken as h/2 
rather than its correct value of (4/2) V1 — (46°/L?); 
for 6/L = 1/10, the error in the result is less than 1 per 
cent. For larger deflections, the influence of this fac- 
tor, together with the fact that the tangent modulus in 
general decreases with increasing load, will eventually 
require a decrease in axial load corresponding to an in- 
crease in deflection. 


BUCKLING OF A COLUMN 


For simplicity consider a rectangular column of 
breadth 6 and depth / (although the method is valid 
Let a set of axes be 
superposed as shown in Fig. 4. At an initial load P,, 
the equation of the centerline is taken to be y = 
y.(x); ata future load P, the equation will be given by 
y = y(x, P). As the load increases from P, to P, let 
the inclination from the horizontal of each cross section 
perpendicular to the centerline change from 6,(x) to 
6(x, P); it is assumed that cross sections initially plane 
and perpendicular to the centerline will remain plane 


for any shape cross section). 


and perpendicular. 
For small deflections, 


d*y, dx? = : Pad = —dé,/dx (14) 


and 
0?y/Ox? = y” = —00/Ox (15) 


As the load increases from P, to P, let the increase in 
compressive strain along the centerline be denoted by 
e(x, P). For those cross sections that do not inter- 
sect any region of strain reversal, e(x, P) will be given 
by 

e(x, P) = (P — P,)/Er bh (16) 


where EF , is the (constant) tangent modulus. 

Consider two cross sections initially Ax apart. At 
load P, their centers will be Ax (1 — €) apart. Points 
on each cross section a distance ¢ to the right of the 
centerline will have their separations decreased by 
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Figure 4 Buckling of a Column 





Sale) = t} A(x + Ax) — 0,(x + Ax) — [0(x) — 0,(x)]} 


The increase in compressive strain at the point (x, ¢) as 
the load increases to P is given by dividing the last ex- 
pression by Ax and taking the limit as Ax tends to zero. 
Using Eqs. (14) and (15) gives 


change in compressive strain = « + f(y” — y,”) 


At any load P, the fibers in the column may be sepa- 
rated into two classes: those whose compressive strain 
has increased steadily during the loading, and those 
whose compressive strain has decreased steadily after 
having increased to a certain maximum. Distinction 
is necessary because the tangent modulus /:; is opera- 
tive for an increase in compressive strain, whereas the 
elastic modulus / must be used for a decrease. In 
order to visualize the stress pattern across a cross sec- 
tion, consider the strain diagram of Fig. 5. The ordi- 
nate represents compressive strain; the distance / 
from the centerline is the abscissa. At some initial 
load, let the strain distribution be AA (since cross 
sections remain plane, the strain distribution must be a 
straight line). As the load increases, the compressive 
strain at the center will increase, and, simultaneously, 
the inclination of the strain line must increase because 
of the larger moment required to balance the effect of the 
larger deflection. It may happen that the inclination 
increases sufficiently rapidly for the compressive strain 
at the outer edge to decrease; such a case is illustrated 
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by BB. The line CC represents the strain distriby. 
tion corresponding to a further increase in load. The 
stress pattern may be obtained by use of the tangent 
modulus for an increase in compressive strain and the 
elastic modulus for a decrease. 

Define the function 1/(t, x, P) by 
ye(x, P) + 


Mt &. PP) = max. 
Po S< Pi <P 
Mh ae >] +. 4 
tly"(x, Pi) — yo"(x)]} (17 


This function represents the maximum compressive 
strain that the fiber at the (x, ¢) position has experienced 
during the loading from P, to P. 
fiber at load P is 


ErM — E{M — [e+ t(y” — y,")]} 


The equations of force and moment equilibrium yi] 


read 


*+h/2 
/ (enM —E{M—[e+t(y’ — ye")]} )bd Is 
h/2 


and 


Py — P.y, = 


°+h/2 
| ies (x rM—E}\ M—([e+t(y"—y,") ]} Yao dt (19 


The stress at this 
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PLASTIC 


If the column is not rectangular, the width 6 will be a 
function of /. Note that the integrands are identical 
save for a factor f and also that the integration is im- 
mediate for any cross section that does not intersect a 
region of strain reversal, since the quantity JJ — 
fe + t(y” — y.”)| then vanishes. Although it is rea- 
sonable to expect the quantity e + f(y” — y,”) at first to 
increase monotonically for increasing P [at fixed (x, ¢) | 
and then to decrease monotonically, this fact has not 
been used in the derivation (for, if stress increases after 
a decrease, the increase is governed by the elastic 
modulus until the previous maximum of stress is 
reached). 

Suppose now that the initial deflection 6, is suffi- 
ciently small that there is no strain reversal for the first 
part of the motion. Then Eqs. (18) and (19) reduce 


at once to 


. | -— ” ” } (20) 
= Feb fe + t(y” — y,”) Jat 
—h/2 


= Faxbh « } 
and 
Py—- Py = — / : Ik ,Mobt dt 
oilitheses *h/2 
= —E,b fe + t(y” — y,”)|t dt) (21) 


—E vb (h®/12) (y” — »,") = 
—E rl (y" _ Vo") 


ll 


Eq. (21) may be rewritten as 


If the initial deflection is assumed sinusoidal, Eq. (22 


iseasily integrated. Let 
Yo = 6, sin (wx/L) (23) 
Then Eq. (22) yields 
y = y,(Pr — P,)/(Pr — P) (24) 


Although this equation was derived on the assumption 
of constant /y, it is valid even if the modulus varies 
during the loading, provided that the deflection is suffi- 
ciently small so that the modulus does not vary appre- 
ciably over the width of the cross section. For if Ey is 


not constant, Eq. (21) may be written 


*h/2 
Py— Py, = — / [Ert(y” — yo") + K|tbdt (25) 

h/2 
where A is a constant equal to the increase in stress at 
the centerline due to the increase in load from P, to P. 
It is assumed in this equation that E, does not vary ap- 
preciably over the width of the cross section. Because 
the constant A cancels out in the integration, Eq. (25) 
reduces again to Eq. (22); since the integration of Eq. 
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(22) is with respect to x, the fact that /; is a function of 
P will not affect the result, Eq. (24). 

Since P 7 decreases with increasing P, the effect of the 
dependence of Ey on P is to make the deflection in- 
crease more rapidly than if /, were constant. Note 
that Eq. (24) is the same as Eq. (5) in the case of the 
Shanley model. 

Expressing the compressive strain at the outer edge 
in terms of P and setting the derivative with respect to 
P equal to zero gives the critical load P, at which strain 
reversal first occurs. If dP, dP is small, the result is 
P, = Pr — V 6(6,/h) P7r(Pr — P.) 26) 


which has the same form as the corresponding Eq. (7) 
in the case of the Shanley model. 

If it is desired to solve Eqs. (18) and (19) for a col- 
umn that has reached the strain reversal load, an ap- 
proximation method must be used. As a first ap- 
proximation to the solution, assume 


e = (P — P,)/Erbh 


and 


y = ¥(Pr — P,)/(Pr —P 


which would represent the solution for no strain re- 
versal. The tangent modulus is assumed constant. 


Let a first approximation to .1/ be given by 

P, — P, 7 (= —_ =) 
; ia : Vo 
E rbh EL? \Pr — P, 


Differentiation of Eq. (27) with respect to P; yields 
that value of P; for which the maximum occurs; the 
first approximation to / is therefore known in terms 
of P, and Py in a strain reversal region. The integra- 
tion of Eqs. (18) and (19) may be broken into two 


M = 


max. 
P «< P; « P 
97 


(26) 


parts: in one part there is no strain reversal and the 
integration is immediate; in the other part the integra- 
tion uses that value of ./ just obtained. The two solu- 
tions must ‘‘line up’’ at the point of demarcation; this 
point is easily found by setting P; equal to P and solv- 
ing for fin the expression previously found, which gives 
that value of P, for which 1/ isa maximum. The final 
result of the integration is that two new values for e and 
y have been obtained; these new values may be used 
as a second approximation, and the whole procedure 
may be repeated. For further details, consult the 
author's thesis. The graphical method sketched in 
the next section is somewhat more direct and incorpo- 
rates the stress-strain curve of the material into the 
calculations [this may also be done with the integral 
Eqs. (18) and (19) but with more difficulty ]. 


GRAPHICAL SOLUTION 


Suppose that at a certain load the deflection of a col 
umn is known at each point, as well as the stress dis 


tributions across certain cross sections. Then, if the 
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stress-strain curve for the material is known, it is pos- 
sible to construct the new deflection curve and the new 
stress distributions corresponding to a small increment 
in load. Repetition of the process will result in ob- 
taining the load-deflection curve of the column. 

It will be again assumed that cross sections initially 
plane and perpendicular to the centerline will remain 
plane and perpendicular. It follows that the compres- 
sive strain must vary linearly with depth across any 
cross section. 

Consider a particular cross section. Let the com- 
pressive strains at the inner and outer extremities of 
this cross section be denoted by e; and «,, respectively, 
and let them be measured from their values at some 
fixed datum load P,. Before the increment in load, let 
the stress distribution across this cross section be repre- 
sented by the curve AA in Fig. 6; here, compressive 
stress is plotted with respect to depth. Let the aver- 
age stress be represented by the horizontal line BB; the 
value of this average stress is, of course, obtained by 
dividing the integral of the stress AA over the cross- 
sectional area by the total area of the cross section. 
If the column were rectangular, the two areas bounded 
by AA and BB would be equal. 

If the load is now increased slightly, let the new 
average stress be represented by CC. Neither of the 
new values of ¢; or €, is known; however, if any value 
for €; is assumed arbitrarily, then the corresponding 
value of ¢, may be found by requiring the average stress 
to coincide with CC. Note that, if the stress curve 
corresponding to the new value of e; intersects the 
curve AA, as shown by DEE, then the actual stress dis- 
tribution must be taken as EEE because of the fact 
that a decrease in stress requires the use of the elastic 
modulus rather than the part of the stress-strain curve 
in question. By a (graphical) integration over the 
cross section, the bending moment corresponding to 
EEE may now be found. Continuing in this way, by 
choosing a sequence of arbitrary values for e;, a rela- 
tionship between e; — €, and bending moment may be 
established for this section. It is necessary to carry 
out this process for a set of evenly spaced cross sec- 
tions—e.g., those marked a, b, c,.... in Fig. 4. 

The relationship between e; — ¢, and curvature 1/p 
must now be investigated. At the datum load P,, 
consider two cross sections a unit distance apart. If 
h is the depth of the column, then the inner and outer 
2p,) and 1 + (h/2p,), where 
If the inner and 


arc lengths will be 1 — (h 
1/p, is the datum load curvature. 
outer compressive strains are now increased by e; and 
€,, respectively, the new radius of curvature is easily 


found to be 
1/p = (1/p.) + [(e, — &)/h] (28) 


Furthur, the bending moment corresponding to any 
value of e; may be equated to the product of column 
load and deflection; thus a relation between e; — €, and 


y is obtained. Use of Eq. (28) then gives the desired 
relationship between 1/p and y at each of the cross gee. 
tions a, b, c, etc. Since 1/p may be approximated by 
y”, the problem is now to solve an equation of the form 
y” = f(x, y), where the functional relationship js, of 
course, known only for certain values of x. The bound. 
ary conditions are y(o) = y(L) = 0 and y’(L/2) = 9, 
Since the solution is symmetrical, only one-half of the 
column need be considered. Several methods for 
solving equations of the above form are given in refer. 
ence 3; a modification of Kelvin’s method (p. 170) will 
be used. 

The deflection of the center is not known, and it js 
therefore necessary to assume some reasonable vyalye 
for this quantity. The radius of curvature may then 
be calculated from the relation between curvature and 
deflection previously discussed. Because the slope of 
the desired solution is zero at cross section a, the center 
of curvature must lie on AA (Fig. 7), say at point A, 
With this point as center, an arc is drawn, intersecting 
the line BB through b at B. The distance Bb is then 
used as the deflection of the column at b; as before, the 
corresponding radius of curvature will be known, 
The first circular are is now discarded, and a new arc 
is drawn with the center of curvature again on AA but 
with the new value for radius of curvature. This are 
is extended twice the distance of the previous arc to 
meet CC at C, and is assumed to represent the deflec- 
tion of the column between a and c. Continuing from 
the point C by the same double jump method, the de- 
flection of the end of the column is finally obtained. 
If this deflection is not sufficiently close to zero (say 
within | per cent of middle deflection from zero), then 
another value of center deflection must be assumed and 
the procedure repeated. The radii of curvature are 
so large that all intersections, etc., must be calculated 
rather than obtained by compass and ruler. This is 
done with the help of Fig. 8, using 


S => (d?/2R cos 6) + d tan 6 


the change in deflection corresponding to any are is 


easily calculated. 


BUCKLING OF A RECTANGULAR PLATE 


The problem of plastic buckling of a rectangular 
plate under compressive edge thrust has been treated 
by Ilyushin,‘ Bijlaard,® Stowell, Kaufman, and Handel- 
man and Prager.’ Allexcept the last of these treatments 
use deformation theories and are therefore open to ob- 
jection.’ The following will be based on the results of 
reference 7; with the addition of a prime, equations 
will be numbered as in reference 7. 

If a rectangular plate of width } and length a is sub- 
jected to a compressive stress ¢, on the two ends, in the 
direction of the x axis, then Handelman and Prager 
show that the criterion for loading or unloading is given 


by 
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PLASTIC BUCELING OF 


W = —(0,/3) (2é — & — é,) = O, respectively 
For loading, their stress-strain relations (in reduced 


form) are 


A-1\. 
Loum. — (tT Cy 


1-1 


= 211 + v) try 


| 


Yr 


where pv is Poisson’s ratio, \ is the ratio of E to Ey, and 
the dot indicates differentiation with respect to time. 
For unloading, the usual elastic relations hold. In 
contrast to reference 7, it will now be assumed that 
lv > 0 over the entire thickness of the plate; this as- 
sumption will be considered in more detail later. Then 
proceeding as in reference 7, the buckling equation be- 


comes 
DuWeese + 2DeWeevy + Dx2W yyyy = —ohWes (72") 
where 
W = deflection rate 
(r 3) (1 — pv’) 
ink Be —— 
(5 — 4v)\ — (1 — 2p)? 
1 (1 — v?) 
Doo = D . ~T 9 
(5 — 4v)X — (1 — 2p)? 
2A — 1 + 2») | 
nO ~ se 
7 ' E — 47)A\—(1—27)*? I+, 
(73") 
D = WE/12(1 — v?) 


h = plate thickness 

It might be objected that Eq. (72’) is not valid for 
use in finding the critical load because the coefficients 
are derived on the tacit assumption that buckling com- 
mences from a state of zero deflection and because for 
such buckling there must be unloading for some por- 
tion of the plate. However, an examination of the 
various assumptions made in deriving Eq. (72’) shows 
that it remains valid for sufficiently small deflections 
and that for such small deflections an initial motion 
of pure loading is possible. For comparison, consider 
the case of a column; here, buckling from a state of 
zero deflection requires a strain reversal (unloading) in 
a portion of the column, and yet the critical tangent 
modulus load may be obtained by finding that load for 
which all deflection rates are possible from an initial 
state of zero deflection, the basic equation being set up 
on the basis of no strain reversal. As in reference 7, 
Eq. (72’) may now be used to find the critical stress 
¢--. for various boundary conditions. 

For a simply supported plate, Eq. (113) of reference 
7 may be written 
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Fig. 10. Comparison with Experimental Results 





Dy = [(\ + 3)/4A]D 
De = Du = D 


and, by differentiation of Eq. (113’) with respect to m 
(the restriction that m be integral is of little importance 
for long plates), the lowest critical stress is found to be 


Cor, = 2(m?D/b*h) [V (A + 3)/4 + 1] 
Since the critical stress for the elastic case is 
Cer = 4(2?D/b2h) 


the ratio 7 = o¢,./¢¢7, will be given by 


n = (1/2) [WA + 3)/4 + 1] 


As in reference 9, define the ‘‘calculated elastic critical 


compressive strain’’ to be 
€cr, = Ger. (1/nE) 


then using the experimental results of reference 9 to 
evaluate \ for various stresses, o,,, may be plotted 
against €,,, as shown in Fig. 10. The experimental re- 
sults of reference 9 are shown as circles. 

Exactly as in reference 7, the limiting case of an ex- 
tremely narrow strip may be investigated; the resultant 
buckling load is found to be 

P... = Eql (x*/a?) (99’) 
which agrees with the tangent modulus theory for the 


buckling of columns. The result does not depend on 


Poisson’s ratio. 


It is of interest to consider the effect of a smal] initial 
deflection, given, say, by 
W =A sin (ry/b) sin (mrx/a) (99 


If the corresponding initial load is denoted by o < g.. 
then the initial equation of equilibrium is 

k 0°?M, oR 07M,, LE 07M, j OV 

q — — 2E q — = oh (3 

Ox? Oxdy Oy" Ox? 7 

where / is the plate thickness and J/,, etc., are the 
usual bending and twisting moments. If the streg 
increases to o + do, let the new deflection be 


W + dW = (A + dA) sin (ry/b) sin (mrx/a) (3) 


The new equation of equilibrium will read 


E OM, D.. dW oR OM, 2D. dW 

~ Ox? sali ~~ Oxdy ee ee em 
O'My : 

E — Dye AWryyy = (6 + do)h(W + dW),, (32 


Oy" 
Subtraction of Eq. (30) from Eq. (32) gives 


—Dyd@Werry — 2DedWery, — DodWyyy = 
h(W,do + odIV,,) (38 


Substituting from Eq. (31), Eq. (83) becomes 


7m \* 7* m? 
dA Di _ 2D 0 a ra ed 
a b? a? 
r\4 7m" 7-m* 
Dey ) — ho — | = hA ——do (34 
b a’ a 


Setting vy = 1/2 as before, Eq. (84) becomes 
dA(o.r, — ¢) = Ada 35 

which may be integrated at once to give 

A = A, (Ger, — 0°) /(Ger, — @) (36 
and this gives the relation between the deflection A and 
the stress o, given that at a certain initial stress o° the 
deflection amplitude was A,. The equation is, of course, 
valid only for small A and A,. 

If it is assumed that unloading occurs over a portion 
of the plate thickness, then, if the load is allowed to in- 
crease during buckling, the methods of reference 7 may 
be used to show that the dividing surface between 


regions of loading and unloading is given by 


l6¢,(v2-1) 


| 
| 
5.4 es | = = 3 (di) 
a + Vv l) eK (4y—S)Eh ~ 


where 


5S,’ = 2/(&/2) 

a =1 — [2/c(5 — 4y)] 

c = f(A — 1)/[(5 — 4x)A — (1 — 20)2]} 
K = (2—») W.- + (2v — 1) Ww 


Now, A ~ A, and Eqs. (35) and (37) show therefore 
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Increased Jet Thrust from Pressure Forces 


F. P. DURHAM* 
Unwersity of Colorado 


SUMMARY 


For ram-jets and turbojets, which operate with low overall jet 
expansion ratios, it is possible to obtain more thrust with incom 
plete expansion in the jet nozzle utilizing both pressure and mo 
mentum forces than with complete expansion utilizing momentum 
forces alone. Proof of this phenomenon is developed, and charts 
are introduced showing the range of pressure ratios and velocity 
coeficients for which incomplete expansion will produce greater 


thrust than complete expansion 
SYMBOLS 


The following notation w il) be used: 


F. = net thrust, lbs 
F, = gross thrust, lbs 
F. = specific gross thrust, Ibs. per lb. air per sec 


= airflow, lbs. per sec 
fuelflow, Ibs. per sec 


i] 


V, = airspeed, ft. per sec 
V, = jet speed, ft. per sec. 
V.' = reversible adiabatic jet speed, ft. per sec. 


= ambient pressure, lbs. per sqft 


p 

p,, = total pressure before jet nozzle, lbs. per sq.ft 

p, = pressure at jet nozzle exit, lbs. per sq.ft. 

1 total temperature before jet nozzle, °R 

T, = temperature at jet nozzle exit, °R. 

7 temperature at jet exit after reversible adiabatic ex 
pansion 

h,, = stagnation enthalpy before jet nozzle, B.t.u. per Ib 

h, = enthalpy at jet nozzle exit, B.t.u. per Ib. 

h,’ = enthalpy at jet exit after reversible adiabatic expansion 

A, = effective jet exit area, sq.ft. 

¢, = specific heat at constant pressure, B.t.u. per lb. °R 

Com = mean specific heat at constant pressure, B.t.u. per Ib 

R 

k = ratio of specific heat at constant pressure to specific 
heat at constant volume 

k,, = ratio of mean specific heat at constant pressure to 
mean specific heat at constant volume 

g = acceleration due to gravity, 32.2 ft. per sec.” 

R = gasconstant, ft.lbs. per lb. °R. (53.3 for air 

J = mechanical equivalent of heat, 778 ft.lbs. per B.t.u 

7; = fluid specific weight at jet nozzle exit, lb. per cu.ft 

@ = V,/V,', nozzle velocity coefficient 

A = correction factor for divergence of nozzle exit 


INTRODUCTION 


[ IS COMMONLY BELIEVED that a jet propulsion power 

plant will produce more thrust with complete ex- 
pansion in the jet nozzle than with incomplete expan- 
sion. However, for turbojet and ram-jet engines, 
Which operate with essentially low-expansion ratios 
across the jet nozzle, it is actually possible to obtain 
more thrust from a converging nozzle utilizing both 
momentum and pressure thrust than from a correctly 

Received January 7, 1950 

* Assistant Professor of Aeronautical Engineering 


designed converging-diverging nozzle that operates 
with complete expansion and utilizes only the mo- 
mentum thrust. The following development is offered 


as proof of this phenomenon. 


DEVELOPMENT OF FORMULAS 


The most commonly used equation for jet thrust is 
obtained from Newton's third law and is 


F, = (wa/g) (V; — Va) (1) 


This equation assumes complete expansion in the jet 
nozzle to atmospheric pressure and considers the weight 
flow of fuel to be negligible compared to the weight 
flow of air. 

Actually, since complete expansion in the jet nozzle 
does not always occur, a pressure thrust term must be 


introduced. This is given by! 


g) + A;(p; — pa) — (WaVa/Z (2) 


This equation assumes that all pressures external to the 
engine except the pressure at the jet exit are atmos- 
pheric. It also neglects the weight flow of the fuel. 

If the weight flow of the fuel is considered, its mo- 
mentum change must be added to Eq. (2), which then 


becomes 


Another factor needed to make the equation com- 
pletely accurate is the nozzle divergence correction 
factor A, which corrects for the fact that the exit veloc- 
ity in a diverging nozzle is not all axial but has a radial 
component.” The introduction of this factor in Eq. 
(3) gives 

bm SUS an (pj) — Pe) -—Ve (4 
£ g 

Since the pressure ratios to be considered are all low, 
a nozzle divergence angle of less than 10° could be used, 
and, consequently, the factor \ can be taken as unity as 
shown in reference 2. In addition, since the introduc 
tion of fuel complicates the thermodynamic calcula- 
tions, especially with regard to specific heats, the weight 
flow of the fuel will be neglected in this study, and the 
development will involve the application of Eq. (2) 
instead of Eq. (3). The principles are perfectly gen- 
eral, however, and can be applied to the flow of com 
bustion gases, as well as air. 

From Eq. (2) the gross thrust F, is defined as 


F = (w,/g) V, +A (p; — Pa) (5) 
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Fic. 1. Effect of velocity coefficient on specific cross thrust. 


F, is a parameter that gives the thrust characteristics 
of the nozzle itself independent of the engine inlet air 
momentum. 

It is more convenient to express Eq. (5) in terms of 
the thrust per pound of air per sec., the specific gross 
thrust, which may be accomplished by dividing Eq. 
(5) by w,. Then 


F, = V; g + (A; Wa) (p; “-~ Pa) (6) 


A further refinement can be made by expressing 
A,/w, in terms of Vj, p;, R, and 7; which may be ac- 
complished as follows: 

From the continuity equation and the equation of 
state for an ideal gas, 


Ww, = A;Vjy; (7) 
1; = p;/RT; (8) 
w,/A,; = RT;/V3p; (9) 


Substituting Eq. (9) in Eq. (6), 
F, = V;/g + (RT,/V;) [1 — (pi/p;)] (10) 


In the case when the jet velocity is sonic, the accepted 
equation is 


V; = a = VkgRT, (11) 


This equation is based on the assumption of reversible 
adiabatic flow.* Since the flow outside the boundary 
layer is essentially frictionless, it approaches a reversible 
adiabatic flow, and Eq. (11) may be used in the actual, 
as well as the ideal, expansion path. 

If Eq. (11) is substituted in Eq. (10) and similar 
terms are collected, 


e = VRTj/kg lk +1 —(fa/p)] — (y 


The subscript ¢ is introduced to indicate a converging 
nozzle expanding to sonic velocity when subjected toa 
pressure ratio greater than the critical. 

The effect of friction on Eq. (12) is shown only in the 
ratio p,/p;, which will be larger in the actual case thay 
with reversible adiabatic expansion. The effect ¢j 
friction in the case of complete expansion is to reduc 
the discharge speed V; 

If a reversible adiabatic expansion is assumed with 
constant specific heat, Eq. (10) may be expressed as, 
function of p; alone for any given overall pressure ratio 
Puj/ Pa, aS Shown in Appendix I. If this expression js 
differentiated with respect to p; and set equal to zerp 
for a maximum, it is found that a maximum exists 
when pj = fp, (that is, when complete expansion 
exists). 

However, in the case of adiabatic flow with friction, 
which exists in real nozzles, the discharge velocity is not 
so large as would be obtained in the reversible adiabatic 
case between the same pressures but can be expressed in 
terms of the ideal reversible adiabatic velocity as fol. 


lows: 
V; = oV;’ (13 


where ¢ is the velocity coefficient and \’,;’ is the revers- 
ible adiabatic jet velocity. 


EXPLANATION OF CHARTS 


Fig. 1 is a plot of gross thrust against velocity co- 
efficient for complete and incomplete expansion for an 
overall expansion pressure ratio, ~,;/p,, of 3.5 and a 
total temperature before the nozzle of 1,600°R. This 
plot indicates that if the velocity coefficient is less than 
0.98, more thrust can be obtained with incomplete ex- 
pansion to sonic velocity than with complete expansion 
The method for obtaining this plot is shown in Appen- 
dix (II). 














Sy 

WAR GREATER || 

S|} RUST WHA +—t 

SS 

e| FT. | 
WJ 

NS 

>| 6 

M 

e| 4h 

= 4 
NS 

§| 71 nore lr} =0p% | | 











/00 098 096 5 OG O90 
VELOCITY COEFFICIENT -Q_ 





Fic. 2.¥ Dividing line between complete and incomplete expat- 


sion for maximum thrust. 
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INCREASED JET THRUST 


Fig. 2 is a plot of expansion pressure ratio against 
velocity coefficient with the points determined from the 
intersections of complete and incomplete expansion 
curves for various pressure ratios plotted in the same 
manner as Fig. 1. To the right of the line in Fig. 2 is 
the region where more thrust would be obtained with 
incomplete expansion, while to the left of the line is the 
region where more thrust could be obtained by com- 


plete expansion. 
DISCUSSION OF RESULTS 


Velocity coefficients of at least 0.95 are attainable.‘ 
With @ = 0.95, it is evident from Fig. 2 that, for overall 
pressures less than approximately 6, more thrust could 
be obtained from incomplete expansion than from com- 
plete expansion. Since this covers the overall expan- 
sion range for turbojets and ram-jets operating under 
Mach 2, design of nozzles with convergence only will 
actually produce more thrust than with converging- 
diverging nozzles designed for complete expansion at a 
given overall pressure ratio. If the velocity coefficients 
of the converging-diverging nozzle is less than 0.95, 
the dividing line in Fig. 2 is not crossed until even a 
higher overall expansion ratio than 6 is attained. 
Thus, if the velocity coefficient is 0.92, more thrust 
will be obtained with incomplete expansion at all over- 


all pressure ratios up to 9. 


It should be noted that the method of Fig. 1 assumes 
the same efficiency for both types of nozzles, while, 
actually, the converging nozzle, if properly designed, 
would probably have a higher efficiency than the con- 
verging-diverging nozzle. Also, if the divergence angle 
is large, there is a further loss due to flow divergence. 
Therefore, the case for incomplete expansion is some- 
what stronger than that shown in Fig. 2. 


As the expansion pressure ratio increases, the differ- 
ence between the ideal thrust with complete expansion 
and the ideal thrust with incomplete expansion becomes 
greater, so that the compensating effect of the velocity 
coefficient is a smaller factor. Hence, above an ex- 
pansion ratio of 6 to 8 with attainable high-velocity 
coefficients, complete expansion will finally produce 
more thrust than incomplete expansion as shown in 
Fig. 2. 

Thus, for the ratios encountered in 
rockets, complete expansion will give more thrust 
than incomplete expansion, with the difference be- 
tween the two increasing as the overall pressure ratio 
increases. This effect for high pressure ratios is well 
known, having been presented by Malina,® Zucrow,® 
and others. While the velocity coefficient is taken as 
unity in these references, the comparison between com- 
plete and incomplete expansion shown is valid, since 
the lowest overall expansion pressure ratio given is 
20, which is well above the favorable range for incom- 
plete expansion as shown in Fig. 2. 
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MAXIMUM THRUST FOR IDEAL EXPANSION 


In the case of reversible adiabatic flow with constant 


specific heat, 


Vy = V2gJ(hy — hy’) = V2gJe(Ty — T;') (14) 
Since 
T'/T 15 = (5/Py)*~“ (15) 
then 
v, i V 2gJe,T yl — (pj/Py)*~» ey 
29 Jc Z. ee net 7 
Toe (py) /k _ (p,)* 1 /k) (16) 
tj 
Now let 
C = V2gI¢T y/(py)*~* (17) 
then 
V,’ a tur" ‘ke g,°-» k) I/s (18) 
Substituting Eq. (19) into Eq. (10), 
C r—]) . (k—1) r\ i/o 
F, = (py ‘_ pM ' ey i/e 4 
RT;{1 — (b./P))] 
7, (k-D7k =D (19) 
Cpr Pi ) 
Now let 
CQ, = T y/(Pu)*~”” (20) 
then, from Eq. (15), 
F, = (py* ko pt» bys 4 
g 
aca? P anu (k-1)/k.-Y/2 mo ; 
ee - (py P= py P“") [1— (be/P)) (21) 
: : RC 
F, = : 1 ities " p;“ —_— = Xx 
‘ j C 
1 alte k _ ites a 2 (p,*-! Jk Dab; ' ky (22) 
and, for maximum F,, 
OF, C | k - 
-=0= +{—-2+21—- babs - 
Op; RCig k — 1 
ee i 
but —D/k _ pt i (9) 
This condition is satisfied when p; = /~,, which cor- 


responds to complete expansion. 
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Appendix I] 
METHOD FOR Fie. 1 


From the general energy equation for any adiabatic 
flow through a nozzle 


Vi = V2gT(hy — hj) (24) 


while, for reversible adiabatic flow between the same 
pressures, the ideal velocity is given by Eq. (14). 
Then, from Eqs. (13), (14), and (24), 


= V (hy — hj)/(hy — h,’) (25) 


J 
In the case of complete expansion for an arbitrary 
overall pressure ratio and an arbitrary nozzle inlet total 
temperature, the reversible adiabatic enthalpy drop, 
h,,; — h,’, may be determined from reference 7. Then, 
for a given value of the velocity coefficient, ¢, the actual 
enthalpy drop, 4,; — h;, may be obtained from Eq. 
(25) and the jet speed from Eq. (24). The specific 
gross thrust may then be calculated from Eq. (10) with 
ba/ ps = 1. This specific gross thrust is then plotted 
against velocity coefficient in Fig. 1 for an overall pres- 
sure ratio of 3.5 and an inlet total temperature of 
1,600°R. 
In the case of incomplete expansion to sonic velocity, 
it follows from Eqs. (11) and (24) that 
2gJ (hij — hj) = kgRT; (26) 
By making use of a mean specific heat 
Com( 23 — 13) = hy — hy (27) 
and the well-known relation 
Jepm/R = km/(Rm — 1) (28) 
where k,, is the ratio of mean specific heats, Eq. (26) 
may be simplified and reduced to 
T15/T; = 1+ (Rm 1) (R/R»)/2 (29) 
However, for expansion to sonic velocity 7, — 7; is 
approximately 300°F., so that the difference between 
k,, for the process and k at the nozzle throat is negligible. 
Substituting k,, = k in Eq. (29) 
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T;, = 2T1/(k + 1) (39 

The value of & for a given nozzle inlet total tempera. 
ture and corresponding pressure ratio (approximate) 
two) for sonic velocity may be readily determined fron 
achart in reference 8. The temperature 7°; can they b 
calculated from Eq. (30). If more accuracy is desired 
new value of k can be obtained from the above cha 
using the inlet total temperature and the value of 7 
just determined. Then 7); can be recalculated frop 
Eq. (30). 


necessary. 


In general, however, this refinement js no 


Using the value of 7’; thus determined, /; may be ob 
tained from reference 7, and the reversible adiabaj; 
enthalpy drop may be calculated from Eq. (25). Thep, 
from this reversible adiabatic enthalpy drop, the cor 
responding pressure drop is found from reference 7 
which permits the determination of p,/ p;. Finally, the 
specific gross thrust can be calculated from Eq. (12) for 
any given arbitrary overall pressure ratio and velocity 
coefficient. These values for an overall pressure ratiy 
of 3.5 and an inlet total temperature of 1,600°R. ax 


plotted in Fig. 1. 
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The Influence of Design Parameters on the 
Performance of Subsonic Air Inlets 


HOWARD E. ROBERTS? anv B. D. LANGTRY? 
Douglas Aircraft Company, Inc. 


SUMMARY 


This paper presents the results of an experimental investiga 
tion into the influence of design parameters on the performance 
of a certain class of subsonic air inlet. This class of inlet, nor- 
mally located in the skin of the air frame, employs a boundary- 
layer bleed. 

The influence of aspect ratio, degree of protuberance, ramp 
profile, boundary layer thickness, and bleed position upon the 
inlet-total-pressure recovery, critical Mach Number, and bound- 
ary-layer-bleed-total-pressure recovery was determined. The 
results are summarized in the form of graphs and design charts 
that can be used by air-frame designers in the selection and de- 
sign of air inlets that satisfy their requirements. 

The application of these results to practical design problems is 


discussed. 
INTRODUCTION 


pag OF modern single-engined military air- 
craft revealed that a typical airplane has more 
than a dozen inlets, and that the number and variety 
of inlets is increasing with the complexity and speed of 
the aircraft. In many cases these inlets must be de- 
signed to suit given requirements, with a reasonable 


assurance of their successful operation, without develop- 
ment testing. Such development testing of the large 


number of inlets of a single airplane is costly. Conse- 
quently, there is a need for the systematic investiga- 
tion of various types of inlets to provide data that can 
The published data of 
The only 


be used as a basis for design. 
such investigations are extremely meager. 
complete systematic design data are for nose inlets.! ° 

The need for air-inlet design information was recog- 
nized by the Power-Plant Division of the Bureau of 
Aeronautics, Department of the Navy, and a contract 
was let with the Douglas Aircraft Company, Inc., for a 
systematic study of air inlets. The present investiga- 
tion was planned to supply design information for sub- 
sonic air inlets, which are normally located in the skin 
of the aircraft and which employ a boundary-layer 
bleed. This paper presents a summary of the results of 
this investigation. 

There is still a critical need for systematic research 
on other types of inlets and for additional research to 


Paper presented at the Aircraft Design Session, Eighteenth 
Annual Meeting, I.A.S., New York, January 23-26, 1950 
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determine the influence of other design parameters on 
the type of subsonic air inlet considered in this paper. 


NOMENCLATURE 


An air inlet is the entrance section of an internal-flow system. 
In this paper an inlet is considered to include the external sur- 
faces in the region of the inlet plane and a short length of inter- 
nal duct. The performance of such an inlet is described in terms 
of its total-pressure recovery and critical Mach Number. This 
performance is influenced by certain design parameters that the 
air-frame designer can control. The symbols listed below are used 
to define the parameters that describe the geometry of the inlet 
and the conditions of the airflow. Some of these symbols are 
shown in Fig. 1 


A = area of duct at total-pressure rake, sq.in. 

a = portion of area, A, sampled by a single total-pressure 
tube, sq.in. 

b = longitudinal distance between bleed-lip leading edge and 
inlet-lip leading edge, in. 

¢ = aconstant used in empirical equation-defining boundary- 
layer profile, dimensionless 

d = duct depth, in. 

H = total pressure, lbs. per sq.ft. 

h = bleed-duct height, in 

ZL = number of total-pressure tubes in rake, dimensionless 

m = an empirical factor used in defining the boundary-layer 
thickness, dimensionless 

n = numerical index, dimensionless 

P = protuberance distance, normal to surface, in 
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Fic. 1. Model configurations and notations. 
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EACH COMBINATION OF GEOMETRIC PARAMETERS INDICATED 
BY SHADING WAS TESTED AT INDICATED COMBINATIONS OF 
BOUNDARY-LAYER PARAMETERS 


Fic. 2. Parameter combinations tested; inlets with circular 
ramp. 


p static pressure, lbs. per sq.ft. 

q = dynamic pressure, lbs. per sq.ft. 

R = ramp radius, in. 

r = inlet-lip leading-edge radius, in. 

t = bleed-lip thickness, in. 

V = flow velocity, ft. per sec. 

X = longitudinal distance from inlet-lip leading edge to point 
of tangency, in. 

= space coordinate, parallel to surface, in. 

Y = distance from inlet-lip leading edge to point of tangency, 
normal to surface, in. 

y = space coordinate, normal to surface, in. 

6 = boundary-layer total effective thickness at inlet location 


in absence of inlet, in. 
x = protuberance function = P/(Ymar. +r +h +t+d) 
(see Fig. 1), dimensionless 


Subscripts 
b = conditions in boundary-layer bleed duct 
n = numerical index 
0 = free-stream conditions 
1 = conditions in main duct a distance d downstream from 
the lip leading edge 
y = conditions at a distance y from the surface 


DESCRIPTION OF EXPERIMENTAL APPARATUS 


The experimental data were obtained from tests con- 
ducted in the 30- by 45-in., low-speed, closed-circuit 
wind tunnel located at the Douglas Aircraft Company’s 
El Segundo Plant. Wooden models of the air inlets 
were mounted in the floor of the test section. An 
auxiliary centrifugal blower was used to draw air 
through the main and boundary-layer-bleed ducts of 
the models. The flow in these ducts was regulated by 
cone throttles and measured by means of flat-plate 
orifices in standard A.S.M.E. test ducts.* Multiple- 
tube manometers were used to measure the model- 
surface static pressures and total-pressure distribu- 


tions. 
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DESCRIPTION OF MODELS AND TECHNIQUE 


Time and cost limitations usually make it impossible 
to include all combinations of geometric and aerody. 
namic parameters in any systematic investigation 
Fig. 2 shows schematically those combinations oj 
parameters that were tested with one of the ramp pro- 
files. In order to investigate as many of the combina. 
tions as possible, the models were designed with inter. 
changeable components so that changes in model cop. 
figuration could be made quickly and simply. 

Tests were conducted to study the effect of the fo. 
lowing geometric variables: degree of protuberance, 
aspect ratio, ramp profile, and boundary-layer-bleed 
position. Fig. 1 shows drawings of the models tested, 
Flush, semiflush, and full-protruding models were used 
to show the influence of degree of protuberance on inlet 
performance. The degree of protuberance is measured 
by the function 7, which varies between zero for a flush 
inlet and unity for a full-protruding inlet. Inlet 
widths of 10, 4, or 2 in. were used to obtain data at 
inlet aspect ratios of 5, 2, and 1, respectively. The 
inlet depth, d, was held constant at 2 in. Previous 
experience indicates that inlets of aspect ratio 5 or 
greater can be considered essentially two-dimensional, 

It will be noted in Fig. 1 that the profile of the outer 
surface of the inlet lip is an NACA Series 1 Nose-Inlet 
Profile.!_ The dimensionless coordinates of such a pro- 
file are given in Table 1 (the lip leading edge is taken as 
the origin of the coordinates). It is recognized that 
a considerable degree of arbitrariness exists in deter- 
mining the relation between the various dimensions of 
the inlet. The relative dimensions of the inlet used 
during the subject investigation are based on typical 
current design practice. Obviously, many other model 
configurations are possible. 


TABLE | 
Dimensions of NACA Series 1 Nose-Inlet.! (Origin of Coordinates 
at Leading Edge) 


/X y/Y a/X y/Y 
0 0 0.400 0.7475 
0.025 0.1657 0.500 0.8269 
0.050 0.2436 0.600 0.8911 
0.100 0.3613 0.700 0.9395 
0.150 0.4530 0.800 0.9735 
0.200 0.5270 0.900 0.9940 
1 


0.300 0.6489 1.000 0000 


Two types of ramp profiles were included in the 
various inlet configurations. The majority of tests 
was conducted using a ramp of circular profile, the 
radius of the circle being 20 in. (R = 10d). A character- 
istic of this ramp is that the angle that the duct axis 
makes with the free-stream direction varies with the 
degree of protuberance. The second ramp used has 4 
profile consisting of a straight section inclined 15° to 
the direction of free-stream flow and of a transition sec- 
tion that is a portion of a circle of 20-in. radius. Inlets 
incorporating this ramp have a constant 15° inclination 
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PERFORMANCE OF 


of the duct axis, regardless of the degree of protuber- 
ance. In each case the duct cross-sectional area was 
maintained constant for a minimum distance of 6 in. 
(3d) downstream of the inlet survey plane. Both 
of these ramps were selected because of their geometric 
simplicity and because they represented practical 
configurations for aircraft applications. 

The boundary-layer-bleed lip was fabricated from 
0.125-in. copper sheet. It was designed so that its 
fore-and-aft location could be varied. Tests were 
conducted with the bleed-lip leading edge located 
0.375d, d, and 1.625d forward of the inlet-lip leading 
edge. The bleed lip was installed in the model so that 
a constant bleed-duct height (i = 0.125d) was main- 
tained. The aft bleed-lip position was not tested with 
the flush inlet, nor the forward position with the full- 
protruding inlet, because the resulting model configura- 
tions were not considered realistic. 

Static-pressure orifices were installed on the ramp, at 
the boundary-layer-bleed lip, and on the inlet-lip inner 
and outer surfaces. These orifices were located on the 
model centerline. Total-pressure rakes were installed 
on the model centerline in the main duct at a distance 
of one inlet depth, d, downstream from the lip leading 
edge (Fig. 1) and in the boundary-layer-bleed duct in 
order to measure the total-pressure distribution. Addi- 
tional static-pressure orifices and total-pressure rakes 
were installed off the centerline to obtain data for de- 
termining the effect of aspect ratio on critical Mach 
Number and recovery. 

The influence of inlet-velocity ratio (Vi/ Vo), bleed- 
velocity ratio (V,/Vo), and boundary-layer thickness 
(6/h) on air-inlet performance was investigated during 
the test program. The boundary-layer thickness was 
determined, with the models removed, at a point on the 
test-section floor corresponding to the leading edge of 
the inlet ramp. This thickness was varied by means 
of roughness elements in the contraction cone of the 
wind tunnel and a wind-tunnel-boundary-layer bleed 
at the upstream end of the test section. The boundary- 
layer thickness is determined by fitting the measured 
velocity profile with an equation of the form 


log y = m (V,/Vo)? + ¢ (1) 


and defining the thickness, 6, as the value of Eq. (1) 
when 


(V,/Vo)* = 1.0, or 
logé6 =m+c (la) 


This may be accomplished graphically by plotting the 
individual rake-tube measurements of (V,/Vo)? vs. 
log y and drawing a straight line through the points to 
(V,/Vo)? = 1. 

The inlet-velocity ratio was varied from zero to in- 
finity on all tests. Tests were conducted at bleed- 
velocity ratios of 0.7 and 0.3. The former value was 
selected because previous tests had indicated this to be 
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near the optimum value. The latter was selected be- 
cause of the smaller boundary-layer-bleed power re- 
The boundary-layer thickness is expressed 


quirement. 
The specific 


as the dimensionless parameter 6/h. 
values of this parameter which were studied were 1, 2, 
and 4. 


REDUCTION OF DATA 


The inlet and boundary-layer-bleed total-pressure 
recovery is defined by the following equations: 


- HH, — Po 
inlet-total-pressure recovery = (2) 
qo 
HH, — po 
bleed-total-pressure recovery = —————_ (2a) 
go 


Eqs. (2) and (2a) are calculated from the data of the 
individual total-pressure tubes according to the follow- 
ing relations: 
. © a, Va (He — po) 
inlet-total-pressure recovery = Set ae eases 
n=1A; Vy Jo 

(3) 
‘25 Ay Ve (H,, ooh Po) 


bleed-total-pressure recovery = : 
n=1A, Vy qo 


(3a) 


Eqs. (3) and (3a) are general and include the effect 
of unequal tube spacing. They are derived by con- 
sidering the flow energy of each fluid element and inte- 
grating with respect to its mass flow. This method of 
computing recovery gives a true indication of the energy 
that is available in the air stream. 

The definition of recovery given in Eqs. (2) and (2a) 
has the disadvantage that, at a value of the inlet- 
velocity ratio of infinity (i.e., when the free-stream 
velocity is zero), the value of the inlet-total-pressure 
recovery is always negative infinity. Therefore an 
inlet-loss factor is used to express inlet-total-pressure 
recovery at extremely high values of the inlet-velocity 


ratio. The factor is defined as 


inlet-loss factor = (JI) — I1)/q (4) 
This definition was selected because the inlet-loss factor 
remains finite when the free-stream velocity is zero. 
It is computed from the total-pressure-tube data by 


the equation 


jue 
a 
i) 
_ 
=. 


inlet-loss factor = — > (fo — Hn) (5) 


The critical Mach Number is determined from the 
static-pressure distribution on the surface of the ramp 
and on the inner and outer surfaces of the main lip. 
The Karman-Tsien method‘ was used in determining 
the values of critical Mach Number presented in Fig. 3. 
Although it is recognized that this method is applicable 
only to two-dimensional flows, it is used here because 
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Fic. 3. Design chart showing inlet pressure recovery and critical 


no adequate theory exists for the type of flow encoun- 
tered in two- and three-dimensional inlets. 
RESULTS 

The results are presented as graphs and design 
charts with the parameters expressed in dimensionless 
form. The main design chart (Fig. 3) demonstrates 
the effect of degree of protuberance and aspect ratio on 
the performance of inlets with circular ramps having 
boundary-layer configurations described by the follow- 
ing values of the boundary-layer parameters: }/d = 
1, V,/Vo = 0.7, and 6/h = 2. Tables 2 and 3 sup- 
plement Fig. 3 and show the effect of the boundary- 
layer parameters b/d, V,/Vo, and 6/4 on recovery and 
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Mach Number as functions of protuberance and aspect ratio. 


critical Mach Number. A comparison of values ob- 
tained from the chart with experimental values for the 
same configuration is presented in Fig. 4. This illustra- 
tion shows that linear interpolation for intermediate 
inlet-velocity ratios yields close agreement with the ex- 
perimental values. 

The inlet-total-pressure recovery is highest for inlet 
velocity ratios near unity. At low inlet-velocity ratios, 
the boundary-layer effects tend to reduce the recovery. 
At high inlet-velocity ratios, there is a contraction of 
the flow, accompanied by internal separation, which 
causes a reduction in recovery. In general, these effects 
are most pronounced for low values of the protuberance 
The influence of ramp 


function, 7—1.e., for flush inlets. 
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Fic. 4. Comparison of design chart data with experimental re- 
sults. Circular ramp, A.R. = 5, = 0.6, b/d = 0.375, Vi/ Vo = 
0.3, 6/h = 4 


profile on the inlet-total-pressure recovery is greatest 
for flush inlets. The straight profile gives better re- 


covery at extremely low values of inlet velocity than 


TABLE 2 
Corrections to Inlet-Total-Pressure Recovery for Effect ¢ 
Boundary-Layer Parameters* 


=~ 


V:/ Vo 0 0.4 1.0 ee: 1.6 
o/h l 0.022 0.020 0.011 0.008 0.008 
2 0 0 0 0 0 

4 —0.004 —0.028 —0.028 —0.026 —0.024 
Vi/Vo 0.3 —(0.166 —0.068 —0.014 —0.012 0 
0.7 0 0 0 0 0 
b/d 0.375 0.008 0.014 0.007 0.007 0.006 
] 0 0 0 0 0 
1.625 0.094 —0.014. —0.016 —0.016 —0.017 


* Corrections to be added (algebraically) to inlet-total-pressure 
recovery determined from Fig. 3 to account for effect of boundary- 
layer parameters. 


TABLE 3 
Corrections to Critical Mach Number for Effect of Boundary- 
Layer Parameters* 


Vi/ Vo 0 0.4 1.0 i. 1.6 
o/h | 0.002 0.002 0.002 0.002 0.002 
2 0 0 0 0 0 

4 —0.002 —0.003 —0.004 —0.004 —0.005 
Vi/Vo 0.3 0.013 0.009 0.004 0.003 0.002 
0.7 0 0 0 0 0 
b/d 0.375 0.021 0.016 0.009 0.009 0.009 
] 0 0 0 0 0 
1.625 —0.007 —0.061 —0.036 —0.028 —0.012 


_*Corrections to be added (algebraically) to critical Mach 
Number determined from Fig. 3 to account for effect of boundary- 
layer parameters 
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does the circular-ramp profile, because of the smaller 
turning of the flow. At all other values of inlet- 
velocity ratio, however, the circular ramp gives the bet- 
terrecovery. (See Fig. 5.) 

The test results show that the recovery is highest 
for the high inlet aspect ratios, because the low-energy 
boundary-layer air along the duct sidewalls is a rela- 
tively small portion of the entire flow. As the inlet 
becomes narrower—i.e., as the inlet aspect ratio be- 
comes smaller—this low-energy air becomes a larger 
portion of the flow, and the recovery falls off. The 
amount of low-energy air along the duct walls is greater 
for flush inlets, because boundary-layer air from the 
tunnel floor can flow over the ramp walls and into the 
inlet, while for the full-protruding inlet the side fair- 
ings prevent this flow. Therefore, the effect of inlet 
aspect ratio on the recovery of full-protruding inlets 
is smaller than for flush inlets. 

The effects of boundary-layer and bleed conditions 
are slightly greater for flush inlets than for the semi- 
flush and full-protruding inlets. As might be expected, 
the recovery is best with the thinnest boundary layer 
and with the greatest bleed-velocity ratio. The re- 
covery is also best, except at low inlet-velocity ratios, 
for those configurations with the aft bleed-lip position. 
At low inlet-velocity ratios, a forward bleed-lip position 
reduces the tendency for the boundary layer to separate 
from the ramp upstream of the bleed-duct lip, and there- 
fore the recovery is improved. 

Inlet loss factors for high inlet-velocity ratios are 
presented in Fig. 6. It should be noted that these 
data are plotted against the reciprocal of the inlet- 
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Fic. 5. Effect of ramp profile on recovery. 








fir AR MM LAM 6 ILA DILE 


sere 


Fl he eHE: 


434 JOURNAL OF THE AERONAUTICAL SCIENCES-—JULY, 1950 






































25 

*]520 

= ¢ AR=5 

a \ b/d= 1 

2 18 

Q VjNet 0.7 

u 6/h=2 

2 ic ” 

°o e " 

- \e a 

© 05 s, va 

a N =95) thd LT. 
—o 

i es cn 

















0 
fe) 2 4 6 8 lO 812° «(14 
RECIPROCAL INLET-VELOCITY RATIO— V9/V, 


Fic. 6. Inlet loss factors for circular-ramp inlets. 


velocity ratio, Vo/V;. The inlet loss is primarily a con- 
traction loss and therefore will be greatest for the 
flush inlet, except at extremely high inlet-velocity 
ratios—i.e., V;/Vo near infinity. The relatively sharp 
leading-edge radius under this condition of flow causes 
the flow to separate. Since this sharp lip extends 
around a smaller percentage of the inlet perimeter on 
flush inlets, the flow is accelerated more gradually, and 
there is a smaller loss. 

The critical-Mach-Number contours in Fig. 3 are 
determined by the inlet component that has the lowest 
critical Mach Number. At different inlet-velocity 
ratios, different components of the inlet have the lowest 
critical Mach Number. For example, at high inlet- 
velocity ratios—i.e., V;/Vo 2 1.2—the inner surface 
of the lip is the component that has the lowest critical 
Mach Number for all inlet configurations tested. At 
these high velocity ratios, the stagnation point is on the 
outer surface of the lip, and, hence, there is a high-veloc- 
ity region around the lip leading edge. At a velocity 
ratio of unity, however, the component that has the 
lowest critical Mach Number varies with degree of 
protuberance. Therefore, the background of the criti- 
cal-Mach-Number contours of Fig. 3 has been shaded, 
to indicate which component has the lowest value of 
critical Mach Number. The ramp is the critical com- 
ponent, at a value of 1’; 1) = 0.8, for wide flush inlets 
(x = 0, A.R. = 5). 
critical Mach Number for the semiflush and _ full- 


The ramp has an extremely high 


protruding inlets. 

At lower inlet-velocity ratios (V1/Vo < 0.8 except as 
noted above), the external lip surface has the lowest 
critical Mach Number. At these velocity ratios, the 
stagnation point is on the inner surface of the lip, and 
again there is a high-velocity flow around the lip lead- 
ing edge. 

The boundary-layer-bleed-total-pressure recovery de- 
pends upon the relative stream energy of the air drawn 
into the bleed duct. In the boundary layer, the stream 
energy varies with the distance normal to the solid 
boundary, while outside the boundary layer the stream 
energy is essentially constant. Therefore, the bleed 


recovery in a duct of constant height will be greateg 
when the boundary layer is thinnest. For a given 
thickness of boundary layer at a point upstream of the 
inlet, the boundary-layer thickness will increase with 
distance downstream and will also increase with th 
degree of curvature of the flow. Hence, it follows tha 
flush inlets, with long ramps and considerable ramp 
curvature, have lower bleed recovery than the ser. 
flush and full-protruding inlets. 
the forward bleed position should give the highest 
bleed recovery, and the aft bleed position should give 
the lowest. Similarly, a thick boundary layer and q 
low bleed-velocity ratio will each cause low bleed re. 


It also follows that 


covery. 
Experimental bleed recovery data are given in Fig, 7, 
These data agree, in general, with the analysis pre- 
sented above. In Fig. 7a, the extremely low bleed- 
total-pressure recovery of the flush inlet at low inlet- 
velocity ratios and the inflection in the recovery curve 
indicate separation of the boundary layer ahead of the 
bleed. In Fig. 7f the recovery curves are in agreement 
with the analysis presented above, except at low inlet- 
velocity ratios where the recovery curves cross, indi- 
cating highest recovery for the thickest boundary 
layer. The reasons for this effect are not known. The 
static-pressure distribution on the ramp _ indicates 
separation at a constant value of pressure coefficient. 
Therefore, the three curves of Fig. 7f should all have the 
same value—equal to the pressure coefficient at the 
separation point—at V,/ Vo = 0, which they do. 


DESIGN APPLICATIONS 


The design charts and illustrations for this paper are 
intended to provide air-frame designers with a rational 
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Fic. 8. 
means for the selection and design of air inlets. Use of 
the design chart for the selection of an inlet requires a 
knowledge of (1) the boundary-layer thickness at the 
proposed location of the inlet, (2) the internal-flow- 
system requirements and pressure-drop characteristics, 
and (3) the required critical Mach Number of the 
inlet. Each of these quantities may be evaluated by 
conventional aerodynamic design analyses. The inlet 
configuration that can best meet the requirements may 
then be selected from the design charts. 

Although this paper is not intended to be a treatise 
on the design of air inlets, a few general comments on 
air-inlet design considerations may be in order. Two 
of the most important considerations in the selection of 
the air inlet are the boundary-layer thickness at the 
location of the air inlet and the required total-pressure 
recovery. For example, in some cases where the boun 
dary layer is thin or the inlet-total-pressure recovery 
requirements are low, the air inlet may not require a 
boundary-layer bleed. Where such an air inlet is 
Suitable, it is preferred to an inlet with a boundary 


UNLESS OTHERWISE SPECIFIED, V./V=l, AR=5, VY/Vo=0.7 


Inlet-total-pressure profiles. 


layer bleed because of its greater simplicity. Avail- 
able space is frequently an important consideration 
in the determination of air-inlet geometry. The circu- 
lar-ramp inlets have the advantage over the straight- 
ramp inlets in that they require a smaller longitudinal 
distance for installation. 

Some types of inlets are required to operate at zero 
forward speed—i.e., infinite inlet-velocity ratio. The 
inlets of a jet engine or auxiliary power unit are ex- 
amples of this requirement. Since the performance of 
each of these power units depends upon the total- 
pressure recovery at the inlet to the unit, the inlet 
loss under this flow condition is important. The data 
presented in Fig. 6 should be useful in estimating the 
inlet losses in the velocity-ratio range corresponding to 
take-off and climb conditions. It is unfortunate that 
the inlet-total-pressure-loss coefficient for the extremely 
high and infinite velocity ratios are rarely presented in 
experimental results of air-inlet investigations. These 
data are extremely important, since they affect the 
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A Theoretical Investigation of Heat 
Transfer in the Laminar Flow 
Regions of Airfoils’ 


LEONARD GOLAND} 
Cornell University 


ABSTRACT 


An analysis is presented of the thermal boundary layer about 
a certain infinite cylinder having a chordwise pressure gradient. 
An incompressible fluid of Prandtl Number equal to unity is con- 
sidered. The thermal boundary-layer profiles are obtained at 
various chordwise stations and are shown to differ appreciably 
from the existing velocity boundary-layer profiles. The method 
of solution is analogous to that used by Sears! in calculating the 
spanwise flow in the boundary layer of a yawed infinite cylinder. 
Although use of this method is not feasible for an arbitrary 
cylinder, it presents a case in which the results of more approxi- 
mate methods may be evaluated. Consequently, the heat 
transfer from the cylinder is calculated by several commonly 
used methods and is compared to the more exact solution. 

The results indicate that, of several available approximate 
methods applicable to the estimation of heat transfer in this 
case, only the one proposed by Squire’ gives good agreement with 


the more exact results. 


INTRODUCTION 


TT PROBLEM OF theoretically calculating the ther- 
mal boundary-layer profiles, about an arbitrary 
airfoil having a chordwise pressure gradient, is at pres- 
ent too complicated for a rigorous solution. There- 
fore, in order to solve the dependent problem of heat 
transfer from an airfoil, various simplifying assump- 
tions have been made as to the shape of the thermal 
layer. Consequently, several approximate methods*~* 
are now available for calculating the heat transfer from 
an arbitrary airfoil. 

It was noted! that the equation for the spanwise ve- 
locity component in a steady, incompressible flow over 
a yawed infinite cylinder is identical to the equation for 
the temperature distribution in a steady two-dimen- 
sional laminar boundary-layer flow, provided that the 
Prandtl Number is unity and the viscous heating is 
neglected. Therefore, the distribution of spanwise 
velocity near a yawed cylinder is the same as the dis- 
tribution of temperature when the cylinder surface is 
maintained at a uniform temperature different from 
the stream temperature. Sears presented a rigorous 
solution to the spanwise velocity gradient distribution 
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about a certain infinite yawed cylinder.' Cong. 
quently, this solution is used here to obtain the therma| 
boundary-layer profiles about the same _ cylinder 
These profiles are then used as a basis to evaluate the 
various simplifying assumptions made in the above 


references. 
SYMBOLS 
Cp specific heat of air at constant pressure 
fe = filmcoefficient = g/(T> — Ts) 
K = thermal conductivity of air 
L- = airfoil chord 
p = pressure at edge of velocity boundary layer 
po = free-stream static pressure 
Py = Prandtl Number = Cpu/K 
gq = rate of heat transfer per unit area 
Re = Reynolds Number = UL/v 
T = temperature at a point in the thermal boundary layer 
Ts = free-stream static temperature 
Ty) = airfoil surface temperature 
R = free-stream velocity 
U = free-stream velocity in x direction 
V = free-stream velocity in y direction 
u = velocity in boundary layer in x direction 


= velocity in boundary layer in y direction 
= velocity in boundary layer in s direction 
velocity at edge of velocity boundary layer in x direc- 


“" = 
tion 

x = distance measured along the airfoil surface from the 
stagnation point in plane normal to the longitudinal 
axis of the cylinder (see Fig. 1) 

y = distance measured along the airfoil surface in the direc- 
tion of the longitudinal axis of the cylinder 

2 = distance measured normal to the surface of the airfoil 
on plane perpendicular to the longitudinal axis of the 
cylinder 

8 = a constant, whose value is dependent upon the thick 
ness ratio of the cylinder 

un = coefficient of viscosity of the fluid 

p = density of the fluid 

v = kinematic viscosity = u/p 

# = anondimensional film coefficient = fc/(AK Re /*/L) 


denotes the value at the surface of the cylinder 


THEORY 
For steady incompressible flow, neglecting viscous 
heating, the equation governing the temperature dis- 


tribution is” 
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The specific heat and thermal conductivity of the fluid 
is assumed constant in formulating the above equa- 
tion. For the case of a yawed or unyawed cylinder in- 
finite in extent in the y direction, Eq. (1) reduces to 


of or 4 w >) = KO 4 —~) 9) 
I an | te Ox? = Oz? “ 
since 0/Oy = O, where x, y, and gz are curvilinear co- 
ordinates defined as in Fig. 1 and wu, v, and w are the 
velocity components in the directions x, y, and 2, respec- 
tively. -By repeating, with the necessary alterations, 
Prandtl’s argument concerning the magnitudes of the 
various terms in the boundary-layer equations, one can 


reduce Eq. (2) to” 
u(07'/dx) + w(OT/dz) = (K/pC,) (077/027) (3) 


In applying Prandtl’s argument, it is necessary to as- 
sume that the radius of curvature of the airfoil is large. 
Therefore, Eq. (3) is not expected to give accurate re- 
sults in the immediate vicinity of the leading edge of 
the airfoil, where the radius of curvature is small. 

It is evident from Eq. (3) that the temperature dis- 
tribution is essentially a two-dimensional problem and 
is independent of the position along the longitudinal 
axis of the cylinder. Therefore, the cylinder may be 
yawed or unyawed, the temperature distribution being 
dependent only on conditions in the xz plane as defined 
herein. 

The spanwise flow boundary-layer equation for a 
yawed infinite cylinder in steady flow is! 


4(Ov/Ox) + w(dv/Oz) = v (07v/0z") (4) 


If we assume a fluid of Prandtl Number equal to unity, 
then K/pC, = v, and the form of Eq. (3) is identical 
to that of Eq. (4). Therefore, a solution of Eq. (4) for 
visa solution of Eq. (3) for 7. Consequently, the dis- 
tribution of spanwise velocity near a yawed cylinder is 
the same as the distribution of temperature when the 
cylinder surface is maintained at a uniform tempera- 
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ture. Sears' considered the case of a symmetrical 
cylinder at zero angle of attack, having an external po- 
tential flow velocity of the form 


uy = BU(E — €') (5) 


A solution to Eq. (4) was obtained for this case and is of 
the form 


>> &G,(n) (6) 


n 0 


where § = x/L andy = (8'R,*/L)z. The boundary 


conditions placed upon G,(m), in compliance with the 
boundary-layer theory, are 

G,(0) = 0, for all n 

G,,( © ) 0, forn = 1 

Gol co) = l 


II 


The values of G,(n) were evaluated and are plotted in 
reference 1. The appropriate solution of Eq. (3) there- 


fore takes the form 


m%—-T=(M—-T) XS &G,(n) (7) 
a ‘= 0 

The shape of the thermal boundary layer determined 
directly from Eq. (7) is shown in Fig. 2 for & = 0.2, 0.5, 
and 0.7. It has often been assumed that the shapes of 
the velocity boundary-layer and the thermal boundary- 
layer profiles about an airfoil are similar. Prandtl* 
has obtained the velocity boundary-layer profiles about 
the cylinder under consideration. His results are also 
plotted in Fig. 2. Note particularly that the slopes and 
thicknesses of the two layers differ considerably. As 
is shown later, this fact explains one of the errors often 
introduced by methods commonly used in solving the 

heat-transfer problem. 
To obtain the rate of heat transfer per unit area, it is 
necessary, in accordance with the classical heat-transfer 
theory, to evaluate the slope of the thermal boundary 


layer at the surface—that is, 
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Thermal and velocity boundary-layer profiles at three 
chordwise stations. 
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Nondimensional heat transfer film coefficient vs. chordwise distance along the cylinder. 
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where the subscript 0 denotes the value at the surface. 


dn 


oT 
(55) 
On 0 dz 


(8) 


From Eq. (7), 


or a - 
-( ) = (To mae i Is & Gy (QO) (9) 
On 0 n 


0 
where the prime denotes differentiation with respect to 


y. Consequently, 


=. eee 
qe kK (170 -—Ts)  &G,'(0) (10) 
4 n= 0 
The film coefficient is defined by the expression 
FAE) = ¢ (To - rs 
so that 
3 : R,' . . 
f(t) = K ~ 2#G,'(0) (11) 
L wpa 
or, letting f.(&) (KR, ”/L) = ® (a nondimensional 


film coefficient), 


y 


n 0 


& = gp” &" G,,'(0) (12) 


For modern airfoils (#/U)maz. iS approximately 


equal to 1.2 at zero lift. For the cylinder under con- 


sideration, 

i} U = pie = &*) 
To obtain a value of (4/U) maz. = 1.2, B is set equal to 
3.0. Taking a value of 8 greater or less than 3.0 would 


increase or decrease (%;/U)mar. from the value of 1.2. 
The potential velocity and pressure distribution for this 


» 


cylinder are shown in Fig. 3. 
» £’"G,,'(0) has been evaluated for certain values 
n 1 

of & and, thus, f.(£) and ®(£) are determined. 

Plot [p] of Fig. 4 shows 

® against £ for the cylinder. It might be noted that 

the separation point of the flow is located at approxi- 

The flow, therefore, is not laminar 


These 


values are shown in Table 1. 


mately & = 0.7. 
past this point, and the above solution is no longer 
valid. 


COMPARISON WITH RESULTS OF OTHER METHODS 


The solutions obtained by use of references 4 through 
‘ are all based upon a fluid having a Prandtl Number 
equal to unity. This conforms with the assumption in 
the foregoing analysis. The results are shown graphi- 
cally in Fig. 4. 


Reference 4 


Pohlhausen found that, for laminar flow along a thin 
flat plate and P, = 1, the velocity and thermal boundary 
layers (as defined herein) are identical.‘ Therefore,’ 


the Blasius solution for the velocity boundary-layer 
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TABLE | 

> é G (Q) po = f 
é n 0 KR, I 
0.1 0.564 0.977 
0.2 0.553 0.958 
0.3 0.522 0.904 
0.4 0.480 0.831 
0.5 0,423 0.733 
0.6 0.344 0.596 
0.7 0.237 0.410 


profile along a flat plate may be utilized in the calcula- 
tion of the heat dissipation. It should be realized, how- 
ever, that this solution is based upon 0p/0x = 0, which, 
of course, is not true for the cylinder under considera- 
tion. 
boundary-layer profile for a flat plate and the existing 


Fig. 5 shows the comparison between the velocity 


thermal boundary layer over the previously analyzed 
0 differ appreciably and, 


cylinder. The slopes at 7 
thus, show the error introduced if the flat plate solution 
is used. From Fig. 5, it is noticed that the slopes (at 
n = 0) of the thermal layers are larger than those of 
the corresponding velocity layers. Therefore, as would 
be expected, the results of this method yield values 
that underestimate the heat dissipation from the cyl- 
inder. This fact is shown in Fig. 4. 


Reference 5 


In reference 5, the airfoil section is considered to be- 
have asa flat plate. Pohlhausen’s solution is used with 
a correction introduced for the existing pressure gra- 
dient. To correct for the pressure distribution, the 
actual fluid velocity near the airfoil surface is used in 
place of the free-stream velocity. However, as is 


shown in Fig. 4, the results obtained by this method 
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also underestimate the heat dissipation from the cyl- 


inder. 


References 6 and 7 


The method of Frick and McCullough® and that of 
Allen and Look’ are similar for a fluid of P, = 1. Funda- 
mentally, these methods also use the Pohlhausen solu- 
tion. However, to account approximately for the effect 
of the pressure gradient existing about the airfoil, the 
thickness of the laminar velocity boundary layer is 
computed by the method of Jacobs and von Doenhoff.’ 
It is interesting to note that, while assuming a Blasius- 
type profile, reference 9 gives excellent results as com- 
pared to the exact results given in references 2 and 3 for 
the thickness of the velocity boundary layer about the 
cylinder. However, the methods of references 6 and 7 
fundamentally assume that the thermal boundary layer 
is similar to the velocity boundary layer. Fig. 2 shows 
that the actual profiles differ appreciably. The results 
obtained by these methods are shown, in Fig. 4, to 
overestimate the heat dissipation from the cylinder. 


Reference 8 

This analysis consists in adopting a standard shape of 
velocity and temperature distribution across the 
boundary layers—namely, the Blasius profile for both. 
The thicknesses of the velocity and thermal layers are 
then determined by use of the momentum and energy 
equations. While this method assumes that the pro- 
files are of Blasius type, an important deviation from 
the other methods previously discussed lies in the fact 
that the thermal layer thickness is considered to be of 
different magnitude than that of the velocity layer. 

The assumption of a Blasius-type profile appears to 
be excellent under all conditions for the thermal layer 
and is good for the velocity layer as long as the pressure 


gradient is small and is favorable—i.e., in the region 
where the flow is laminar. The agreement between 
this method and the method presented is excellent, ex. 
cept for the region aft of the minimum pressure point of 
the airfoil (¢ = 0.576). As noted above, agreement 
cannot be expected in this region. However, since the 
flow after the minimum pressure point is not likely to be 
laminar, this is not a serious criticism. 

Although this method involves a more lengthy com. 
putation than is required for the other methods, it pre. 
sents the most rational analysis to date for the approxi- 
mate solution of the general problem. 
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Note on Free Turbulent Flows, with Special 
Reference to the Two-Dimensional Wake’ 


G. K. BATCHELOR" 
Trinity College, Cambridge, England 


SUMMARY 


For some time past there has been a growing realization that 
the mixing-length theory of turbulent shear flow is neither em 
pirically nor logically sound. However, without considerable 
guidance from experimental research it is difficult to replace it. 
This paper reviews some of the measurements concerning free 
turbulence which have been made by Townsend at Cambridge in 
recent years and which exhibit several weaknesses of the mixing- 
length theory. It is possible to infer from the measurements the 
general mechanism of transfer of momentum and of quantities 
like turbulent energy and heat, although a corresponding analyti 


cal theory has not yet been formulated. 


(I) INTRODUCTION 


_— MIXING-LENGTH THEORY of turbulent transfer 
has had a long and useful career, but it is now 
generally believed to have serious limitations and in- 
its usefulness is slightly 


consistencies. Even 


spurious; the successful predictions made on the basis of 


past 


the mixing-length theory have been concerned with 
the distribution of some simple mean quantity, usually 
velocity, and there now exists evidence that a number 
of equally successful predictions of mean velocity can 
be made on the basis of different assumptions about the 
When an attempt is made to pre- 
say, the mean energy 


turbulent transfer. ! 
dict more than the mean velocity 
of the turbulence--the mixing-length theory definitely 
fails. 

These difficulties show that it is important that we 
have a clear understanding of mechanical processes 
occurring in a shearing motion. The mixing-length 
theory is an analytical formulation of an assumed 
physical picture. Recent experimental work shows 
that this physical picture cannot be correct, so that there 
is no point in seeking different analytical formulations 
of this same physical picture. The need, therefore, is 
first to understand how the transfer occurs and then, 
and only then, to express the transfer analytically so that 
the equations can be solved to give any desired quan- 
tity. 

Experiments have not yet gone far enough to provide 
a complete understanding of the transfer process, but 
a good beginning has been made. This informal note 
is intended to outline some of the current ideas about 
free turbulence (i.e., a turbulent shear flow not influ- 
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enced by rigid boundaries). These ideas and the ex 
periments on which they are based are the work of 
Townsend** at the Cavendish Laboratory. Most of 
the experiments have been made with a two-dimen- 
sional wake, but it is believed that the mechanics is 
essentially the same in all cases of free turbulence. 
The mixing-length theory is examined first in order to 
determine the source of its inadequacy. 

(II) MIXING-LENGTH THEORY 


POSTULATES OF THE 


The ideas stated by Prandtl® are too well known to 
need description. The mixing-length theory of turbu 
lent transfer can be reduced essentially to the following 
three distinct postulates (appropriate to a field of mean 
velocity U’ which varies in the y-direction only) : 

Rate of transfer of any quantity @ across unit 
—e(d6 dy), where ¢€ is a 


(a) 
area normal to y-direction = 
turbulent transfer coefficient. 

(b) «¢ « /v, where / and v are a length and velocity 
characteristic of the turbulent motions or eddies. 

(c) v « lidU/dy|, and soe « /?\dU/dy). 

Postulate (a) is equivalent to the assumption that 
the transfer is carried out by fluctuating motions that 
are small in length compared with the lengths represent- 
ative of the distribution of 6—1.e., that the transfer is a 
“gradient” type of diffusion. Postulate (b) is, of course, 
uncontentious until further the 
magnitude and distribution of / and v are made, and 
postulate (c) makes such an hypothesis. Prandtl® 
based postulate (c) on the assumption that the mean 
velocity is a transferable quantity in the sense of postu- 
late (a) and that the ratio of #2? and v? is approximately 
An alternative interpreta- 


hypotheses about 


independent of position. 
tion of postulate (c) is that it assumes that the turbu- 
lence has a similar structure at all positions and that 
diffusion of turbulent energy is negligible. For the rate 
of conversion of mean flow energy to turbulence 
energy [— uv(dU/dy)] is then equal to the dissipation 
of turbulent energy, which is proportional to (u?) ‘*// 
on the assumption of similarity, and postulate (c) 
follows, since uz, “2, and v? are all proportional because 
of the similarity. 

The theory is not useful until a speculation about the 
dependence of the mixing-length / on position is made, 
but, since this speculation of necessity varies with the 
kind of flow under discussion, it will not be considered 
Does 


here. What does concern us is the question: 
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Fic. 1. Distribution of turbulence intensity across wake. 


there exist a choice of / (consistent with the general 
ideas of the theory) which leads to accurate results? 


(III) DEFECTS OF THE MIXING-LENGTH THEORY 
Some of the important defects are: 


(1) Inreality € is not zero when dU/dy = 0, and the 
theory leads to profiles that are too pointed in the neigh- 
borhood of dU/dy = 0. This is, of course, one of the 
penalties for neglecting [in postulate (c) above] the 
spreading of turbulent energy by diffusion. Prandtl 
recognized the defect and suggested® the use of 


« = [2[(dU/dy)? + 1/°(d2?U/dy?)?]'” 


where /’ is another undetermined length. However, 
the lack of intuitive information about /’ makes this 
modified expresston of doubtful value. 

(2) As remarked above, postulate (a) implies that / 
is small compared with the dimensions of the mean 
flow. A check on the consistency of the theory is ob- 
tained by using it to deduce / from measurements of 


U. In this way it is found that 


1 = 0.66y [where U(y) = (1/2)U(0)] for a two- 
dimensional wake 
= 0.22y [where U(y) = (1/2)U(0)| for a two-di- 


mensional jet 
= 0.4y near a plane wall 
= (.13a at the center of a circular pipe 


In none of these cases is the consistency check satisfied 
adequately. This reveals the inaccuracy of the basic 
assumption of a gradient type of diffusion—i.e., one in 
which / is small. 

(3) The theory takes no explicit account of the pe- 
culiar edge conditions (i.e., the juxtaposition of regions 
of fully developed turbulence and of undisturbed flow) 
which must occur in all cases of free turbulence. In 
the absence of experimental data, this is perhaps not 
an unreasonable assumption. However, the data now 
available (see Section V) show that it is not possible to 
ignore the edge conditions. 

(4) The theory leads to inaccurate predictions of the 
distribution of turbulent energy. The single curve that 


2, and w? varies as (dU/dy)* and has 


it predicts for u?, 2 
the general shape shown in Fig. 1. 


9 


tributions of u?, v?, and w? are shown for compari- 


The measured dis- 


son. 


(IV) LATER MopIFICATIONS 


In 1942 Prandtl’ and Gértler* advocated the replace. 
ment of postulates (b) and (c) above by 


€< b( U mes. = a) 


in cases of free turbulence, where 6 is the total width oj 
the mean motion and (Umar. U min.) is the total vari. 
ation of mean velocity across the section of width } 
The idea behind this postulate is that the big eddies ar 
responsible for transfer in free turbulent flows and that 
these big eddies will have a representative length and 
velocity proportional to 6 and (Umar. — Unmin.). We 
note that this idea is inconsistent with the retention of 
postulate (a), which specifically assumes that the trans 
fer is carried out by small eddies and that the transfer 
is proportional to the local gradient d6/dy. Gortler 
claims that better agreement with measured distriby 
tions of mean velocity is obtained with this modifica 
tion of the theory. 

A further modification, without restriction to cases 
of free turbulence, was put forward by Prandtl’ in 1945 
He proposes to replace postulate (c) by an equation for 
the balance of turbulent energy, the quantity (uw? + 92+ 
w*)’* then being used in place of v in postulate (b). 
In principle, this is a great improvement on the orig- 
inal theory, because it eliminates the crude approxima- 
tion that the local turbulent energy is determined wholly 
by the local work done by the Reynolds stresses. How- 
ever, the energy balance equation can be written down 
only when assumptions about all three relevant me- 
chanical effects (dissipation, work done by Reynolds 
stresses, energy diffusion) are made, and it cannot 
be said that the final result is free from arbitrariness. 
As yet, this modification of the theory has only been 
applied to turbulent pressure flow between parallel 
planes, for which there is insufficient experimental 
data to provide an adequate comparison with theorett- 
cal predictions. 

Both of these modifications retain the original as- 
sumption that the transfer is proportional to d6/dy [pos- 
tulate (a)], and we have seen that this assumption 
leads to internal inconsistencies concerning the magni- 
tude of /. This postulate, that the small eddies are 
responsible for most of the turbulent transfer, can be 
shown experimentally to be untrue in one case at least 
viz., the turbulent wake of a long cylinder. 


(V) GENERAL INFERENCES FROM MEASUREMENTS IN 
A Two-DIMENSIONAL WAKE 


The experiments Townsend has made on the turbu- 
lent wake of a long cylinder provide useful information 
about the mechanics of the motion. A considerable 
amount of the research has already been described i 
published papers,’ so that it is necessary only to sum 
marize the general inferences from the experiments. 

(1) These and other measurements show®* that the 
small-scale components of the turbulence have the kind 
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of similarity predicted by Kolmogoroff. That is to 
say, the small-scale components of the turbulence are 
isotropic and have statistical properties that are deter- 
mined wholly by the two parameters, v(viscosity) and e 
(energy dissipation per unit mass). A useful conse- 
quence of small-scale isotropy is that the expression for 
energy dissipation reduces to 15v(Ou/Ox)* (provided the 
Reynolds Number is high enough), so that the dissipa- 
tion can be determined from a single hot-wire measure- 
ment. 

(2) As already remarked in Section III, the quan- 
tities v2, v2, and w® (x and u are parallel to the main 
flow and z and w are parallel to the span of the cylinder) 
are unequal in general. Hence, a satisfactory theory 
must discriminate between the 
It seems reasonable to suppose that v? is diffused later- 
ally more rapidly than vu? and w? (1.e., that v? >> vu?, 
ow?), since high values of v are accompanied simul- 
taneously by high values of the transporting velocity 
This supposition is 


three components. 


and of the transported quantity. 
consistent with the qualitative features of the dif- 
ference between observed distributions of v? and of u? 
and w?; the distribution of v? has fewer peaks and hol- 
lows and extends further into the undisturbed stream. 

(3) Theoretical discussions of motions in which 
the mean velocity is only approximately unidirectional 
(asin wakes and jets) usually assume that the distribu- 
tion of mean quantities across the flow is similar at 
different distances downstream. The experiments 
show that, in a two-dimensional wake behind a circular 
cylinder, this assumption is valid only at great distances 
(beyond about 600 diameters) from the cylinder.” The 
lack of similarity is not easily recognized from the meas- 
urements of mean velocity but is made clear by meas- 
urements of the turbulent energy. The ratio of the 
energy of the mean motion to the energy of the turbu- 
lence increases with distance downstream and slowly 
approaches the value 4.5. Apparently, the disturb- 
ance produced by the cylinder is such as to create a 
greater proportion of turbulent energy than is required 
for the ultimate balance with the mean flow. Other 
cases of free turbulence may not resemble the wake 
closely in this respect. 

(4) Itis easily seen from oscillograph records of the 
velocity fluctuation near the edge of the wake that the 
fluctuation is intermittent—i.e., for a certain fraction 
of the time the velocity varies in the ordinary irregular 
manner, while for the remainder of the time the velocity 
fluctuations are slow and of small magnitude. A 
quantitative description of the intermittency can be ob- 
tained by postulating a sharp division of the flow field 
into a laminar flow ‘outside a wholly turbulent wake 
core with an irregular boundary within which the turbu- 
lence has small-scale isotropy and similarity. Then, 
at any point in the wake, laminar and turbulent flow 
will occur alternately as the irregularities of the core 
boundary are carried downstream. If turbulent flow ‘ 
occurs, on the average, for a fraction y of the total 
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time, measured value of (Omu/Ox)*/[(Ou/Ox)?|? = 


(1/ y) measured value in wholly turbulent flow = 3.5/7. 
Measurements of (Ou/Ox)* [(Ou/Ox)*]? at 
points in the wake thus determine y, and a number of 
similar determinations have been found to give a con- 
More recently,°® 


various 


sistent variation of y across the wake.* 
direct determinations of y have shown agreement with 
the distribution inferred in the above manner. Fig. 4 
shows a typical variation of y compared with the dis- 
tribution of mean velocity. The point to be noticed 
is that the intermittency is not a phenomenon that can 
be ignored, as is implicitly assumed in the mixing 
length theory, and has to be taken into account even at 
positions as close to the wake center as the position of 
maximum shear. 

(5) When allowance is made for the effect of inter 
mittency on the measured mean values, it appears that 
the turbulent core of the wake is approximately homo- 
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geneous.‘ For example, (uw? + v? + w?)/y is approxi- 
mately constant across the wake; likewise, the effective 
momentum transfer coefficient uv/y(0U/Oy) is approxi- 
mately constant. Evidently, the core is composed 
chiefly of eddies that are small compared with the core 
width, and the energy associated with these small 
eddies is spread uniformly over a sharply defined re- 
gion. There is evidence that, on the other hand, the 
bulges at the edge of the core have dimensions com- 
parable with the core width. The turbulence appears 
therefore to be roughly divisible into small eddies that 
comprise the core and large eddies that distort the 
boundaries of the core. The distinction between the 
two ranges of eddy sizes is sharper at higher Reynolds 
Nuinbers. 


(VI) THE MECHANISM OF TRANSFER 


On the basis of the physical picture described briefly 
in the preceding section, it is possible to make some 
tentative inferences about the mechanism of transfer 
across the wake.* > Further work on this problem is 
still in progress, so that the statements made hereunder 
should be treated as sub judice. 

In the first place, we can definitely see that a gradient 
diffusion type of transfer (as envisaged in the mixing- 
length theory) is not the only type of transfer which 
operates. A comparison of measured distributions of 
u? and uv across the wake (as in Fig. 5) shows that 
there is a (small) region where u2v and Ou?/Oy have the 
same sign, so that there the transfer of u? is up the gra- 
dient. Whatever assumptions may be made about 
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the mixing-length /, Prandtl’s theory always leads to , 
positive transfer coefficient and, hence, to a transfer 
directed down the gradient. Similar comparisons oj 
measured distributions of v* and v*, vw? and w?, vT and 
T (T denotes temperature) have been made, and in each 
case it would be difficult, if not impossible, to accoun 
for the relation between the two curves on the basis of 4 
gradient type of transfer. The inference is that, jy 
addition to a gradient mechanism, some other trans. 
porting process is present. 


Now consider the picture of the wake as composed 
of a core whose interior is filled uniformly with smajj 
eddies that contain most of the energy of the turby. 
lence and the boundaries of which are distorted by large. 
scale eddies. In general, the effect of the large eddies 
is to move the core boundaries outward (giving the in. 
crease in width of the wake with distance downstream 
This provides a powerful mechanism of convective 
transfer for quantities that tend to be strongly asso- 
ciated with the core, such as turbulent energy and tem- 
perature (if the cylinder is heated). If the average 
outward velocity of the bulges at the boundary of the 
core is V, the rate of transfer of a quantity 6 associated 
with the core will be simply V6. Again it is found that 
this mechanism of transfer, associated with the large 
eddies in the wake, does not account for the whole of 
the measured transfer. It seems that the small eddy 
or gradient type of transfer and the large eddy or con- 
vective type of transfer make important contributions 
to the transfer of quantities like turbulent energy and 
temperature which are associated with the core of the 
wake. Thus, the total measured transfer of some prop- 
erty 6 at a point in the wake may be written® 


v0 = V0 — K,(00/dy) 


where A, is the eddy diffusion coeflicient describing the 
effect of the small eddies. Recent work® has been 
directed toward establishing numerical values of V and 
K, for various quantities specified by @. As is to be 
expected, A, does not have the same value, on the one 
hand, for quantities like turbulent energy which are 
not dissociated from the transferring motions and, on 
the other hand, for quantities like temperature which are 
dissociated. 


The picture is different in the important case of the 
transfer of mean flow momentum. Momentum can be 
transmitted by pressure forces, as well as by direct ex- 
change of fluid particles, and we do not expect meat 
flow momentum to be associated strongly with the 
wake core. A check that the mean velocity does not 
vary sharply at the edge of the core is provided by 
measurements of 


~ 


L= / f(r)dr 


where f(r) is the longitudinal correlation coefficient 


of the downstream velocity fluctuation. If the meat 
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velocity did vary sharply at the core boundary, there 
would be large-scale fluctuations in the velocity re- 
corded at a fixed point near the edge of the wake, and 
the measured value of L would be large there, whereas 
in reality L does not vary much across the whole 
wake. 

Hence, in this case, the outward movement of the 
wake core produced by the large eddies does not make 
an important contribution to the transfer, and the small 
eddies thus dominate the transfer of momentum. As- 
suming that the effective momentum transfer coefficient 
within the core (y) is constant across the wake (and 
also constant downstream in regions of similarity), the 
effective coefficient at a fixed point is yn, and so we 


should have 
uv = —yn(OU/Oy) (U = mean velocity) 


The measured values of uv, y, and LU do, in fact, satisfy 
this equation with a value of n that is constant across 
the wake. Moreover, if the value of n so determined is 
written as the product of a mixing velocity comparable 
with V v? and a mixing length, the latter is found to be 
small compared with the wake width, which is consist- 
ent with the supposition that the small eddies deter- 
mine the value of 7. 

It seems, then, that the mixing-length theory of 
Prandtl is in error, inasmuch as it does not take proper 
account of the effect of the large eddies. In the case of 
the transfer of momentum, the error is indirect. The 
small eddies of the turbulence are responsible for most 
of the transfer of momentum and this transfer is indeed 
describable by a diffusion coefficient type of theory, but 
the transfer occurs only within the core—i.e., within 
boundaries that are determined by the big eddies. 
Fig. 6 shows the improvement in agreement between 
calculated and measured mean velocity distributions 
when this big eddy effect, or intermittency, is taken into 
account. The three lots of calculated points are based 
on (1) Prandtl’s old mixing-length theory (giving a 
profile that is too pointed at the center) (+), (2) the 
revised theory in which the momentum transfer coeflici- 
ent is assumed to be constant across the wake (®), and 
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Fic. 6. Mean velocity in wake (far downstream ). 


(3) a momentum transfer coefficient 
to the measured intermittency y (4). 

In the case of quantities that are strongly associated 
with the core (turbulent energy, temperature, etc.), 
there are both direct and indirect errors in the mixing- 
length theory. Contrary to the assumptions of this 
theory, both the large- and the small-scale components 
of the motion make important direct contributions to 
the transfer. Also, in assessing the transfer produced 
by the small eddies, it is necessary to take account of 
the fraction y of the total time for which these small 
eddies are in action at a fixed point, as in the transfer 


proportional 


of momentum. 
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Determination of Frequency Characteristics 
from Response to Arbitrary Input’ 


E. R. WALTERS? anp J. B. REAT 
Douglas Aircraft Company, Ine. 


SUMMARY 


A method is given for determining the frequency response 
spectrum of a simple dynamic system when an arbitrary driving 
force (input) and the corresponding response (output) are known 
The complex frequency response (transfer function) is determined 
as the ratio of corresponding terms in the Fourier expansions of 
the driving force and response. The numerical evaluation of the 
Fourier coefficient integrals is shown to be equivalent toa simple 
matrix multiplication. The actual computation may be done 
rapidly on automatic calculating equipment, and an example 


has been carried out. 


INTRODUCTION 


ANY DYNAMIC SYSTEMS, particularly those that 
M are composed of a number of simpler subsystems, 
may be studied conveniently by means of the transfer 
function concept. For example, in the analysis and the 
synthesis of a servomechanism, it is desirable to be able 
to determine its output response to any expected in- 
put. The output should be stable, and the degree of 
stability and the response time should be satisfactory. 
It is, of course, impractical to analyze the system for a 
large number of inputs, so, as a compromise, analyses 
for special or particular inputs are made. Consistent 
with the impedance theory that has been developed ex- 
tensively in connection with the design of electrical 
circuits, the most common operational inputs are a 
sinusoidal function and a step function. The response 
to a sine function is particularly convenient to calcu- 
late, and the response to any input can readily be cal- 
culated if a sine function response spectrum is avail- 
able. This leads to the concept of the transfer func- 
tion, sometimes known as the vector performance 
operator, which relates the output to this particular 
input. A considerable amount of work has been done 
recently with the application of transfer function meth- 
ods to aeronautical problems. References 1 and 2 
describe some of the work done at Cornell Aeronautical 
Laboratory, Inc., and at the Instrumentation Labora- 
tory at Massachusetts Institute of Technology; these 
Received August 11, 1949. Revised and received January 2, 
1950. 

* The authors wish to thank L. G. Campbell, Jr., for carrying 
out the example and J. R. Lowe for the description of the IBM 
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papers have extensive reference lists, and the first has 
a synopsis of dynamic research. 

These transfer functions (performance operators 
are particularly convenient to use in the analysis and 
synthesis of a servo system, using methods such as the 
Nyquist criteria to determine the stability of the sys. 
tem. They do not, however, yield information that 
is as easy to visualize as the transient response, and for 
this purpose the actual time response to a selected in 
put is needed. Generally, the step function or square 
wave is used as a reference input to determine impor 
tant characteristics of the output, such as the degree of 
stability and the time for the transient to damp to some 
given fraction of the steady-state value. This time is 
usually referred to as the response time and is often 
taken as that necessary for the input to damp to, and 
remain within, 5 per cent of its steady-state value. It 
is convenient to have a method for calculating the 
transient (time) response from steady-state (frequency 
response data, and vice versa. It is the purpose of this 
paper to present a method that has been found to be 
convenient and practical. 


(GENERAL APPROACH 


There are two aspects to the problem: Given the 
steady-state frequency response data (transfer func- 
tions), find the transient response to a given arbitrary 
input; or, given the input and the corresponding out- 
put, find the transfer functions. In the first case, it 1s 
only necessary to represent the input as a Fourier 
series, multiply the amplitude and shift the phase of 
each component as indicated by the transfer function 
for that frequency, and sum the resulting components. 
In the second case, it is possible to express both the in- 
put and the output as Fourier series with the same fun- 
damental frequency, and the ratio of the terms of cor- 
responding frequency is then the transfer function at 
that frequency. A nonperiodic function cannot be ex- 
panded in Fourier series directly, but any finite part ol 
itcan. The Fourier integral can be used for some non- 
periodic functions but must be modified so that it be 
comes the Laplace transform for many functions of 
practical interest. The latter involve improper (com- 
plex) integrals and so are less convenient than the 
simple series, particularly when the function to be 
expanded is known only empirically. In the interest 
of simplicity, the series is used here, expanding only 
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DETERMINATION OF 


that part of the function which is affected by the trans 


jent. 

In most of the previous studies, these Fourier series 
have been written for a particular input, usually a 
square OF triangular (pulse) wave, and other forms 
have been approximated by sequences of these par- 
ticular waves. However, the series can be obtained 
directly, though in most cases the integration to ob- 
tain the Fourier coefficients will have to be done by 
numerical or mechanical methods. A tabulating ma- 
chine method based on matrix multiplication has been 
found to be convenient for this numerical integration 
and will be explained in more detail after a brief dis- 


cussion of the basic theory. 


THEORY 


Any ordinary periodic function can be represented by 
a Fourier series. If the independent variable is ¢ and 
the function is f(t), then the exponential form of the 


series 1S 
fi a os > ca ( 1 ) 
where 
l iad P inwot 9 9 
aed f(tye “i (n=, «1, *2,... (2) 
je a : 
and w) = 24/7 (7 is the period of the function f(¢) 


and f is the initial value of f). 

The exponential form of the series is used here for con- 
venience in the following development of the transfer 
functions. Nonperiodic functions can be handled by 
use of the Fourier integral or the Laplace transform, 
but it will be sufficient for the present purpose to con- 
sider only a part of the nonperiodic function. If the 
part is long enough for the function to approach closely 
its steady value, it will contain all the pertinent infor- 
mation. This useful part of the function can be ex- 
panded in a Fourier series in a number of ways; the 
following is convenient and has a simple physical inter- 
pretation. Take the fundamental period to be 7, such 
that the transient has damped to a satisfactorily small 
amplitude at a time less than (1/2)7. To complete the 
cycle, invert the first part and shift it vertically and to 
the right so that it joins the end of the first partial cycle. 
The process is illustrated in Fig. 1 and means physically 
that the moving element is imagined to return to its 
original position after the transient has disappeared 
and in the manner in which it was originally displaced. 
Eqs. (1) and (2) again give the required series. If the 
period 7 is divided into two equal parts, only half of the 
coefficients will be nonzero, because of the symmetry. 
To get both odd and even harmonics, the two parts of 
the cycle must be taken of different duration. Points 
between the original ones may be obtained by changing 
the fundamental period. 
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DEFINITIONS 


T(t = input to system 

O(t = output from system 

R(w = amplitude ratio at angular frequency w 
o(w = phase shift at angular frequency 

wo = 27/7) = fundamental angular frequency 

F(w) = Re’? = transfer function 


When the input to a device is given, it can be repre 
sented as in Eq. (1) by 


n + « 


I(t) = D> J(n)e'"™™ 3) 
where 
Y we 
J(n) = | I(t)e~'"*" dt (4) 
T Jo 


fora = 0, #1, =2..... 


The transfer function may be written R(w)e'®’* 
Then the output component for the mth frequency will 


be 


J(nye™™" R(w)e'*” = J(n)R(n)e\o' re 


since this accomplishes the required amplitude multi 
plication and phase shift. The frequency w in ques 
tion is the mth harmonic, so the transfer function is 
Now 


summing over all input frequencies as in Eq. (3), the 


identified by the coefficient 1, as well as by w. 


total output is 


n + 


O(t) = » I(n)R(n)e™ or re 5 


n 
where the /(m) is determined by the input as per Eq. 
(4). 

The procedure for the inverse situation is as follows: 
Specifically, let the output from a device and the input 
to it be given; then determine the overall frequency 
response (i.e., transfer) function. 

The output may be written as a series, just as the in- 


put was in Eqs. (3) and (4): 


n + x 
O(f) = pe P ser 6 


 & 
P(n) = a O(f)e iment yy a 


land 7 is the same period as in Eq. (4)]. 


where 


The output is, also that given by Eq. (5). Thus, 
n — ” 4 
Of) = >» [J(n)R(n)e? n Je" Pics ys P n ein g 


n 


Now 7, the period, is the same in both Eqs. (5) and 


(7), so, of course, the fundamental frequency wo is the 
same in both. 
Since Eq. (8) must hold for all values of the variable, 


the two series must be equal term by term. Thus, for 


any M, 
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J(n)R(n)e’?” = P(n) 


and referring to Eqs. (4) and (7), 


T 
inwot 
P(n) J ve eeee dt 
= : (9) 


Fie) = Rine™'” = = i 
J(n) — 1mwol 
I(tye "iat 
0 


The transfer function at the frequency mw is thus 
determined for any nm. The F(w) for a frequency be- 
tween two integral multiples of the originally selected 
fundamental wp can be calculated, if it is required, by 
using a slightly different wo. It should be noted that 
the period 7 must be the same in both expansions, since 
Eq. (9) depends on this condition. The relative ac- 
curacy of the calculated transfer function for a par- 
ticular frequency will depend on the proportion of 
this frequency in the given input and output. 

The matrix method for numerically evaluating the 
integrals of the required Fourier coefficients is described 


in the next section. 


MATRIX METHOD OF CALCULATING THE COEFFICIENTS 


This method is readily adaptable to the use of auto- 
matic calculating machines and the timesaving is worth 
while even if only a few sets are to be calculated. The 
trigonometric series is more suitable for computation, 
since the terms are periodic. The relation between the 
coefficients in the exponentials and trigonometric forms 


1S 
Cc, = (1/2) (a, — 40,) (10) 


The coefficient of the mth cosine term 1s 


‘2’. 
an = a | S(t) cos norot dt 


Dividing the period into 1/ + 1 equal time intervals, 
Ai, the coefficient of the ath cosine term is, approxi- 
mately, 


= M 


At m : 
a, = > t(mAt) cos (nwom At) 


= = (10a) 
7 m=0 


and, similarly, for the coefficient of the mth sine term, 


\/ m=M 
b, = Tr > f(m Ai) sin (rwom At) 


(10b) 


Except for the factor At/7, Eq. (10a) can be inter- 
preted as the mth element in a column matrix [a,,], which 
is obtained by the matrix multiplication 


igen! Lf m | _ [a,] (lla) 


(n= 0, 1, 2,3,....N) (m= 0, 1, 2,... Mf), whereg,,, = 
cos (nmMw At) is the element in the mth row and mth 
column of the rectangular matrix [g,,,| and where 


tm = f(mAt) (the mth element of the column matrix) 


is the ordinate f(t) at ¢ = mAt. Similarly, 
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[kam] Lim] = [bn] (1th 


where /,,,,, = sin (nmwo At). 

The index m runs from 0 to M and the index » j 
the number of the coefficient. The number of col 
umns in [g,,] and [h,,,] and the number of rows (el 
ments) in [f,,| are the number of points in the ny 
merical integration. Since the index n gives the har 
monic number of the coefficient, it must run from 0; 
the highest harmonic required. 
tained in the integration only when 1 is considerab) 
The two matrices can be combine) 


Good accuracy js ob 


smaller than m. 
into one, and it is of some convenience to have it }y 
square, as the form (12a) below will be if 2V = Mand, 


runs from 0 to NV. 


nm Qn 
f n = ( Do 
i A | | he a 


Eqs. (lla) and (11b) ean also be combined in the fol 
lowing way, which gives the coefficients of the com- 


plex series corresponding to Eq. (2): 


[gum — thnam| lfm] = [Gn tby| = 2[en}) (12 


Defining the matrix, 


A - foe — tien, | 


(This is the A used in the description of the IBM cal 
culation for the example.) Then, Eqs. (4) and (7) be 


come, respectively, 


J(r) = (At/2T)A[/(m)] | 
P(n) = (At/2T)A[O(m)] § 


(12 


A numerical example is given in the next section to 
demonstrate the application of this method. 


EXAMPLE 


An input (Fig. 1) and a transfer function (Fig. 2a or 
Figs. 2b and 2c) were chosen arbitrarily, and the cor 
responding output was calculated and plotted in Fig 
3. This provided the normal starting point for the 
problem in which the transfer function is desired, the 
input and output functions being known, and allowed 
a comparison of the calculated results with the original 
transfer function. 

Calculation of the transfer function starts with the 
tabulation of a set of points representing the known in 
put and output curves (Tables 1 and 2). The tabula 
tion of the input curve is consistent with the replot 
shown in Fig. Ib, although in actual practice it 1s 
The period of the 
wave in Fig. lb is 8 sec., and it is divided into SO equal 
parts (0.l-sec. intervals). The ordinates of the input 
and output curves at the division points are then tabu 
lated in order. 

The coefficients of the Fourier Series that describe 
the input and output curves can then be calculated, 
using IBM equipment to carry out the two matrix 
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SOLID LINE INDICATES ASSUMED ARBITRARY 
TRANSFER FUNCTION, © POINTS CALCULATED 
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FIG. 2b TRANSFER FUNCTION PHASE ANGLE 


SOLID LINE INDICATES ASSUMED TRANSFER 
FUNCTION, © INDICATES POINTS CALCULATED 
FROM INPUT AND OUTPUT. 


SOLID LINE INDICATES ASSUMED TRANSFER 
FUNCTION, © INDICATES POINTS CALCULATED 


FROM INPUT AND OUTPUT. 
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TABLE 2 
Tabulation of Output Ordinates 
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TABLE | 
Tabulation of Input Ordinates 

Time Time 

(Sec.) Ordinate (Sec. ) Ordinate 
0 0 4.1 0.675 
0.1 0.325 4.2 0.250 
0.2 0.750 4.3 —(.100 
0.3 1.100 4.4 —0.2385 
0.4 1.235 4.5 —(0.248 
0.5 1.248 4.6 —0.200 
0.6 1.200 4.7 —0.065 
0.7 1.065 4.8 0.042 
0.8 0.962 4.9 0.057 
0.9 0.942 5.0 0.061 
1.0 0.940 5.1 0.048 
1.1 0.953 5.2 0.019 
L.2 0.983 5.3 —0.009 
1.3 1.008 §.4 —0.015 
1.4 1.015 5.5 —0.012 
1.5 1.012 5.6 —0.008 
1.6 1.007 a —0.003 
| 1.003 5.8 0.000 
1.8 1.000 5.9 0.002 
1.9 0.998 6.0 0.001 
2.0 0.999 6.1 0.000 
2.1 1.000 G.2 0.000 
2.2 1.000 6.3 0.000 
2.3 1.000 6.4 0.000 
2.4 1.000 6.5 0.000 
2.5 1.000 6.6 0.000 
2.6 1.000 6.7 0.000 
2.7 1.000 6.8 0.000 
2.8 1.000 6.9 0.000 
2.9 1.000 7.0 0.000 
3.0 1.000 0.000 
3.1 1.000 oe 0.000 
3.2 1.000 4.3 0.000 
3.3 1.000 7.4 0.000 
3.4 1.000 7.5 0.000 
3.5 1.000 7.6 0.000 
3.6 1.000 ce i 0.000 
ref 1.000 7.8 0.000 
3.8 1.000 7.9 0.000 
3.9 1.000 8.0 0.000 
4.0 1.000 


multiplications of Eq. (12c). 
multiplying AB = C is now described briefly. 


A machine process for 


The matrix A is punched in a ‘reference’ or master 
deck in four fields as follows: field a, row (n); field 
b, column (m); field c, cos [nm(2/40)]; field d, sin 
[um(r/40)|. It is kept in sequence by 0. 

The column matrix of ordinates, B = [f(mAt)], is 
key-punched from a list as follows: field 6, row (n); 
field e, ordinates [ f(t) |. 

Then B is match-merged in front of A by 6. The file 
is reproduced, reproducing a, b, c, and d and ganging 
e. The new cards are sorted to a, and blanks are 
merged behind them by a. The file is multiplied on 
the 604, accumulating =ce and Sde, which are punched 
on the blanks; a is ganged to the blanks in the same 
operation, and the elements of C are tabulated. If 
a row of check sums is appended to A and handled just 
as any other row, then the sum of the elements of C 
must approach zero, providing a good check on the 
work. 

Having the coefficients /() of the input series and 
P(n) of the output series, the transfer function can be 
calculated by Eq. (9). These coefficients have been 
listed for each multiple frequency in Tables 3 and 4. 
If a, and b, are the real and imaginary parts, respec- 





Time Time 

(Sec.) Ordinate (Sec. ) Ordinate 
0 0.000 4.1 0.909 
0.1 0.141 4.2 0.702 
0.2 0.348 1.3 0.329 
0.3 0.721 4.4 —(). 133 
0.4 1.183 4.5 —(). 503 
0.5 1.553 4.6 —().618 
0.6 1.668 4.7 —(). 459 
0.7 1.509 4.8 ~).110 
0.8 1.160 4.9 0.213 
0.9 0.837 5.0 0.375 
1.0 0.675 5.1 0.339 
La 0.711 5.2 0.16] 
2 0.889 5.3 —(). 082 
1.3 1.182 5.4 —().173 
1.4 1.223 5.5 —(). 193 
1.5 1.243 5.6 —(). 086 
1.6 1.136 &.7 0.075 
7 0.975 5.8 0.112 
1.8 0.9388 5.9 0.102 
1.9 0.948 6.0 0.082 
2.0 0.968 6.1 0.030 
2.1 1.020 6.2 —(). 025 
2.2 1.075 6.3 —(0).045 
ye 1.095 6.4 —().014 
2.4 1.064 6.5 0.022 
2.5 1.028 6.6 0.050 
2.6 1.000 6.7 0.050 
2.7 1.000 6.8 0.012 
2.8 1.038 6.9 —().013 
2.9 1.0638 7.0 —(.023 
3.0 1.073 ° ie 0.000 
3.1 1.050 7.2 0.000 
3.2 1.050 7.3 0.000 
3.3 1.050 7.4 0.000 
3.4 1.050 40 0.000 
3.5 1.050 7.6 0.000 
3.6 1.050 ‘oe 0.000 
3.7 1.050 7.8 0.000 
3.8 1.050 i S 0.000 
3.9 1.050 8.0 0.000 
1.0 1.050 

TABLE 3 
Input Fourier Coefficients 
Frequency 
b i = nfo 

n_ (Cycles per Sec.) n bn 

l 0.125 — 1.55421 25. 66852 
3 0.375 — 1.88605 9.06128 
5 0.625 —2.74175 5.80733 
7 0.875 —3.84747 

9 1.125 —3.56179 

11 1.375 — 2.11080 

13 1.625 —1.14671 

15 1.875 —().71330 

17 2.135 —0.51465 

19 2.375 —(). 45695 

21 2.625 —().39281 

23 2.875 —(). 18521 

25 3.125 —().03739 

a7 3.345 —().04324 —(0. 00666 

29 3.625 —().11683 0.03940 

31 3.875 —(). 15422 0.01900 

33 4.125 —(). 15424 —(.01001 

35 +.375 —(). 12742 —().01633 

37 $625 —0.12173 0.00270 

39 1.875 —(). 13307 0.004538 

All even coefficients are zero except ado = 40.0070 


tively, of the input component for the mth multiple o! 
the fundamental frequency (w = mw) and if p, and 
g» are the real and imaginary parts, respectively, of the 
corresponding component of the output, then the trans- 
fer function at the angular frequency w is 
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FIG.3 OUTPUT DETERMINED BY ASSUMED 
INPUT AND TRANSFER FUNCTION. 


oa P(n) Pn — 14n = 
T(w) = = "= R(w)e’? 
Jin) a, — 1b, 


R(w) = V pr? + Qn?/V a,” + 5,? 
Pw) = tan : (—n Pn) — tan , (—6, Qn) 


The input and output were made periodic with the 
return part of the motion beginning at (1/2)7,, so, be- 
cause of this symmetry, the even terms do not appear. 
The calculations for various odd multiple values of a 
are shown in Table 5, and the computed points of the 
transfer function are plotted in Figs. 2a, 2b, and 2c. 
It will be noticed that the points depart from the orig- 
inal assumed transfer function (solid lines in Fig. 2) 
in the vicinity of nm = 13 (1.6 cycles per sec. approxi- 
mately). However, this departure is not serious 
because the components at frequencies higher than 
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TABLE 4 
Output Fourier Coefficients 


Frequency 





f = nfo 

n (Cycles per Sec.) Pn Qn 

1 0.125 —2.72441 26 .41785 
3 0.375 —3.20885 9.90284 
5 0.625 — 4.58830 6.09430 
7 0.875 — 7.36990 2.61342 
9 1.125 — 5.87002 — 5.13296 
11 1.375 2.35649 —3.22034 
15 1.875 0.02109 0.13897 
17 2.125 0.26571 0.25587 
19 2.375 —0. 14808 0.05203 
21 2.625 0.02599 0.02025 
23 2.875 —0.23341 0.07530 
25 3.125 0.09155 0.00110 
27 3.375 —(0.16614 0.18892 
29 3.625 —(). 14398 —(.12699 
31 3.875 0.03613 0.02016 
33 4.125 —(). 19287 0.07446 
35 4.375 —(.02818 —0.03555 
37 4.625 — (0.00229 0.02364 
39 4.875 —(. 152338 0.03738 

All even coefficients are zero except po = 42.00000. 


original input and output. It may also be noted that 
the points established by this one calculation are too 
far apart to ascertain the exact height of the peak in the 
amplitude ratio plot (Fig. 2c), so some judgment must 
be exercised in plotting. Additional points could be 
obtained by using a different period. 


DISCUSSION 


The present theory is primarily concerned with sys- 
tems that are linear in the time derivatives, since a 
transfer function for more complicated systems would 
not be a function of frequency alone. In many cases, 
it is reasonable to assume that the system is linear over 
the range of variables of practical interest. In fact, 
systems are often designed to be as nearly linear as 
possible for simplicity of analysis, when this is possible 
and satisfactory. In other cases, the above process 
may give a useful first approximation when the non- 
linearity is not too great. 

The harmonic content of the input is a natural limi- 
tation on the accuracy of the method. At frequencies 
where the input has only small components, the out- 
put will also have only small components for normal 


transfer function magnitudes. The calculated transfer 


|.6 cycles per sec. make only small contributions in the function is the ratio of these small components whose 
TABLE 5 

° -Transfer Function 

Frequency Phase 

(Cycles per Sec.) —Output Coefficient -Input Coefficient Amplitude Angle, 

n f = nfo Pn — Qn an —bn Ratio deg 

0 0 1.05 7 0.0 1.00 7 0.0 1.05 0 
l 0.125 —2.724 — 126.418 — 1.554 — 725.669 1.08 —2.5 
3 0.375 —3.209 —7 9.903 —1.886 —7 9.061 1.12 —6.6 
5 0.625 —4.588 —7 6.094 —2.742 —i 5.807 1.19 —11.7 
7 0.875 —7.370 —7i 2.6138 —3.847 —7z 3.380 1.52 —21.7 
9 1.125 —5.870 +2 5.132 —3.562 —7 0.756 2.14 — 53.2 
1] 1.375 2.356 +12 3.220 —2.111 +72 0.445 1.85 —114.3 
13 1.625 1.032 +7 0.134 : —1.147+i 0.479 0.84 —149.9 
15 1.875 0.021 —7 0.139 —0.7138 +172 0.344 0.18 — 236.4 
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accuracy is limited by the accuracy of the given input 
and output functions. 


This method has been found to give good results with 
entirely reasonable computation time. For example, 
10 coefficients can be fitted to 8O points with the tabu- 
lating equipment in 2 hours, which is about the time 
required to determine one coefficient using a desk cal- 
culator. This estimate assumes that the basic matrix 
[Znm — thnm| = A is available as a file of cards, since 
it can be used to expand any function provided the same 
number of points is used. To use a fixed A matrix 


(i.e., a fixed number of points) efficiently, the period 7 


should be no longer than necessary to ensure that the 
transient has been dissipated. 


It is interesting to observe that the sum given by 
Eq. (10) is the same as that required if the expansion 
were made by numerically evaluating the Duhamel in- 


tegral of a step function expansion. It is also the Same 


as required by a process developed by R. C. Seaman 


in which the function is represented by a sequence oj 


triangular pulses. 
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Note on ‘‘Interference Between Wing and Body 
at Supersonic Speeds’”’ 


Ray £. Bolz 
Assistant Professor of Aeronautical Engineering, Rensselaer 
Polytechnic Institute, Troy, N.Y 


Acril 17, 1950 


: er NOTEIS TO CORRECT the second equation of Eq. (7) in the 
article “Interference Between Wing and Body at Super- 
Theory and Numerical Application,’’ by Carlo 


sonic Speeds 
AERONAUTICAL 


Ferrari appearing in the JOURNAL OF THE 
Scrences, Vol. 15, No. 6, p. 317, June, 1948 


If Eq. (5 
10 m 70 
$,,°” = ( ) / f(x ar cosh u) du 
ror -/ cosh x/ ar 
of the above-mentioned article is expanded for m = 2, then 


“0 
1071 : : 
p = f(x — arcosh u) (—a cosh u) du 
ror - cosh v/ar 
os 0 

iffa 

= f(x — ar cosh u) cosh u du + 
rir? 

7 cosh x/ar 


7) 
f(x — ar cosh u) cosh? u in| 


r 
cosh ©! x/ar 


Integrating the first term by parts and combining with the second 


term gives 
cosh x/ar- 


*() 
rp) = @? / fo(x — ar cosh u) (2 cosh? u — 1)du 
° 


This would be the form of the second equation of Eq. (7) in the 


reference article. The second equation of Eq. (8) then becomes 


*0 
0 . 2 ; 
& | = —aq’ fo(x — aR cosh u) X 
or r=R cosh v/ ar 


. 


(2 cosh? x — cosh u)du 


ind the corresponding second equation of Eq. (8’) must then read 


te x12 ; 2x)? x? ; ie (* 
abky! — 4. “~ a ah 
3 \Na?R;? aR,2\ aR? R 


ind, finally, the second equation of Eq. (9) becomes 


(= —x ) ; 
a ae 
a xj —1\? (* —- oe ) —— oe 
7 alt. \ alt, 
Xn, — xi'\? Xn — x;’\2 FSi 
a — —-1l1=-F, ) 
o aR, yv( aR, ) (# 


Correspondingly then the third equations in Eqs. (8’) and (9) 


are, respectively, 


) ( ” ) x,* 
a 3 — 
aR,/ \ a2R,? 


= in => Xj os Xn — x; 2 
Bel VOSEy—- 

l aR \ aR, 

aR, V\ aR, 


The change in these equations involves no error in the Ferrari 
numerical analysis because his numerical example used only the 
¢:'* and k,\” terms and did not involve the second equations at 
all. However, numerical application of the theory using the 
second equations must employ the corrected forms. 
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On the Stability of the Boundary Layer Over a 
Cone 


R. H. Battin and C. C. Lin 

Department of Mathematics, Massachusetts Institute of Technology, 
Cambridge, Mass 

April 19, 1950 


| se THERE HAS BEEN considerable interest in the sta 
bility of the boundary layer over a cone in supersonic flow. 
Since Hantzsche and Wendt! have shown that a suitable trans- 
formation of the equations for the steady flow in the boundary 
layer over a cone leads to those for a flat plate, it is generally 
recognized that the work of Lees and Lin for the two-dimensional 
case? should be applicable to the cone problem as well. How- 
ever, a detailed analysis has so far been lacking. 

The equations for the small disturbances in the boundary layer 
over a body of revolution have now been examined, and it has 
been found that the equations are indeed identical with those for 
a two-dimensional boundary layer, within the limitation of usual 
approximations. Thus, the local stability characteristics of the 
boundary layer over a body of revolution can be examined by 
existing theories. It is only to be remembered that comparison 
must be made for the same Reynolds Number based on the thick- 
ness of the boundary layer. Indeed, the other dimensionless 
quantities must also be defined with the boundary-layer thick 
ness as the reference length. 

In applying the theory to the downstream development of a dis 
turbance of a given time frequency, the case of the cone is ob 
viously different fromthe case of a flat plate, because the bound 
ary-layer thickness varies with downstream distance in different 
manners. Since comparison must be made at corresponding 
Reynolds Numbers based on the thickness of the boundary layer, 
it must involve different distances downstream. In fact, if 
x; and x2 are two stations along the flat plate and if 7, and re are 
two stations along the cone having corresponding Reynolds 
Numbers, the result of Hantzsche and Wendt shows that 


to — % = 3(X2 — X)) 
Based on this result and the present analysis, it can be shown 
that, if the amplification of a disturbance between two stations 
x, and x over a flat plate is A, then the corresponding amplifica 
tion (between 7, and rz), in the case of a cone, is A*. In other 
words, a plot of log A versus x for the flat plate and the cone 
differ from each other only by scale factors of 3 for both coordi 


nate axes. 
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Complementary Energy Method in Vibration 
Analysis 


A. |. van de Vooren and J. H. Greidanus 
National Aeronautical Research Institute, Amsterdam, Holland 
March 21, 1950 


ly THE May, 1949, 1ssue of the JOURNAL OF THE AERONAUTICAL 
ScrENcES (Readers’ Forum, p. 316) E. Reissner outlined a 
procedure—called ‘‘complementary energy method’’—for cal- 
culating frequencies and modes in vibration problems. This 
method was expected to yield more accurate approximations for 
the frequencies than the conventional Rayleigh-Ritz procedure if 
in both cases the same assumed modes are employed. This state- 
ment was recently checked and confirmed by P. A. Libby and 
R. C. Sauer (JOURNAL OF THE AERONAUTICAL SCIENCES, p. 700, 
November, 1949) for the symmetric vibrations of a free-free beam 
with a concentrated mass at its center. 

It may be of interest to remark that the complementary energy 
method (abbreviated as c.e. method) is essentially identical to a 
method established by Grammel! in 1939. Grammel proved 
that, for self-adjoint systems, this method always yields results 
for the frequencies that lie between the exact values and the 
slightly greater values obtained from a Rayleigh-Ritz (R.R.) 
analysis with the same assumed modes. 

Restricting the investigations to bending vibrations for sake of 
simplicity, the set of equations determining the frequencies w 

p 
with corresponding modes Yo aivi(x) is given by [compare Eq. 
j=! 
(5b) of Libby and Sauer’s note | 


— f LM;,(x) M;,(x) 
> a;- w? : d 
, 0 EI (x) 


eT : 
‘— f p(x)¢;(x)exn(x)dx 7 = 0, 
0 f 
eo soe ) 


where \/; represents the bending moment due to the load pg;, 


el, iL wg A 
M;(x) = f dn p(é)e\(E)\dE -| (f: — x)p(£:) 9; (Ei) dé, 
x n x 


(2a) 


1.¢e 


and, similarly, 
PL 
M;,(x) -| (fo — x) p(k) px (E)dbs (2b) 
x 


Substitution of Eqs. (2a) and (2b) into the first term of Eq. (1) 
1 1 


leads to 


< * Lie » » ys * » . » 
; i (1 — x) p(E1) 9) (E1) dé, } ; S. (& — x) p(E) gu (Es) dk} he 
0 EI (x) 
(3) 
The integration is seen to extend over the volume of the pyra- 


mid OABCD of Fig. 1. 
x as the first instead of the last operation, Eq. (3) may be re- 


Taking the integration with respect to 


placed by 


L L <fi, 2 
(é — x) (& — x) 
J “i : dxt x 
0 0 \ 0 EI(x) f 


p(E1)p (Ea) 9; (Er) ee(Eo)dkide (4) 
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where So = ® indicates that the boundaries of integration ar 
are 
0 and the smaller of the two values & and £». 
It is immediately seen that 


<&, aes — x) (& — x) 5 
: dx = G(é, &) 
J 0 EI (x) 
represents the influence function for a bending beam. The set 
of Eqs. (1) then can be written in the form 


2 ( _£L FL : 
> a; 27 FJ, eHleeo(G)G (Es, Eo(Eei(Edede — 
i=] 
L ) 
(x)p(x)¢;(x)dxe = 0,k = 1,2 . 
a G K(x) p(x) pj (x dx ¢ ,k ES a 


It will be shown that these equations can also be obtained in 
another way. When the iterated functions 


PL 
g,%(x) = f, G(x, &)p(E) ej (E)dé 6 


are taken as assumed modes and substituted in the differeptia) 
equation, the result is 


p 2 2.(1)\ 
d? dp; 
p>: aj ; (ar a ) — w*p(x) 9") (x) ; = ¢€(x) 7 


p 
where e(x) is the error arising because >> ajg;\) will, in gen 
j=! 
eral, admit no exact solution. The error, weighted in accordance 
with the original modes ¢;(x), is now made to vanish. Hence, 


Pi, 
Pe g(x) e(x) dx = O 8 


Substitution of Eq. (7) into Eq. (8) leads to Eqs. (5), as can easily 
be checked. 

The last method of deducing Eqs. (5) shows that they differ 
from the Rayleigh-Ritz equations only by taking ¢;‘!)(x) instead 
of ¢;(x) as assumed modes, keeping the weight functions ¢g,(x) in 
Eq. (8) unchanged. It is well known that the transformation 
(6) reduces the components of modes with higher frequencies 
compared with the component of the fundamental. Since the 
components of the p + 1th and higher modes are responsible for 
the errors in the result, it will be obvious? that the c.e. method 
yields more accurate values for the frequencies than the R.R 
method with the same assumed modes. 

To attain actually the improved accuracy for ws’, ..., w,°, it 
will often be necessary to perform the calculations to a greater 
number of figures, since otherwise even desired modes may disap 
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pear by the transformation (6). In general, the number of 
figures in ¢; 1) must be increased with logio w,?/a:? compared with 
Still more accurate results may be obtained by applying the 
RR. method (assumed modes and weight functions identical) 
). The pertaining equations have the 


with assumed modes ¢ 


form 
p Pz 
, 4 dy? | on (x) p(x) ej) (x)dx — 
a a J0 
= nm 
en? (x)p(x) ¢)(x)dx ¢ =(Q (9) 
0 


avoiding again all differential operations. It can be shown that 
in first approximation the errors resulting from the c.e. method 
are geometrically intermediate between those of R.R. analyses 
with assumed modes ¢, and ¢;", respectively. 

Comparing the amount of work when applying the c.e. method 
and the R.R. method with functions ¢;‘, two cases must be dis- 
tinguished. 

(1) The elastic properties are given as a stiffness distribution 
El(x). Application of the c.e. method in the form of Eqs. (1) 
requires for the numerical evaluation of the first term pn + (1/2) X 


p(p + 1) operations of the type 


eh 
f r(x)s(x)t(x)dx 


a 


when » represents the number of stations at the beam. For 


evaluation of the first term in Eqs. (9) 
(1/2)n(n + 1) + pn + (1/2)p(p + 1) 


of these operations are necessary, which usually will be consider- 


ably more. 

(2) The elastic properties are given as influence functions. 
Then the c.e. method in the form of Eqs. (5) requires the same 
number pn + (1/2)p(p + 1) of operations as the R.R. method, 
Eqs. (9). Hence, in this case, the latter method should be pre- 
ferred because of its greater accuracy. 

The foregoing considerations are also valid for systems with 
concentrated masses, which can be included in the mass distribu- 


tion p(x) by means of the singular Dirac function 
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Determination of Flutter Frequencies 
and Speeds 


E. Reissner 

Department of Mathematics, Massachusetts Institute of Technology, 
Cambridge, Mass. 

May 5, 1950 


— VAN DE VOOREN AND GREIDANUS are correct in stating 
that the procedure outlined by me! for the determination 
of flutter frequencies and speeds becomes equivalent to a method 
proposed earlier by R. Grammel when specialized to vibration 
problems in vacuum. I regret not having been aware of this 
fact (which recently was also brought orally to my attention 
by Dr. H. A. Wood, of Chance Vought Aircraft) when writing 
my note on the subject. 

It may be of interest to point out that the approach in my 
note was such as to bring out the connection between the pro- 
cedure in question and the well-known minimum principle for 
the stresses in static elasticity, which, it seems, had not been 
done previously. It was, in fact, the quest for an extension of 
known results of this nature which led to my note. In this 
connection, it may be appropriate to make reference to another 
note on the problem of vibrations of three-dimensional elastic 


bodies.? 
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Bifurcation Criterion and Plastic Buckling of Plates and Columns 


(Continued from page 424) 


that, if A is sufficiently small, no solution can exist for 
S,"—1.e., loading must take place over the entire thick- 
ness of the plate. The deflection at which the solution 
breaks down may be obtained from Eq. (37); beyond 
this point, numerical methods are necessary. 
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The Influence of Design Parameters on the Performance 
of Subsonic Air Inlets 


(Continued from page 435) 


take-off distance of jet aircraft by their influence on 
the engine thrust. 

Air-frame-compartment vents are a special type of 
air inlet. These vents operate at low or zero inlet 
velocities and are required to maintain a certain pres- 
sure. It can be seen in Fig. 3 that the total-pressure 
recovery available at zero inlet-velocity ratio depends 
primarily upon the degree of protuberance. The 
higher the required inlet recovery, the greater the 
protuberance required. 

Boundary-layer-bleed-duct design frequently pre- 
sents a number of problems. Since the bleed-total- 
pressure recovery is comparatively small, it is ideally 
suited to ventilation of compartments or cooling of 
accessories, which require a small pressure drop. Wher- 
ever possible, this air should be exhausted through a 
jet ejector or in the region of the jet nozzle in order to 
reduce the internal drag of the boundary-layer-bleed 
system. If this is not possible, however, the boundary- 
layer-bleed air may be exhausted through carefully 
designed exits in the surface of the air frame. 

A study of the drag of air inlets is beyond the scope of 
the present investigation. The determination of air- 
inlet drag and its subdivision into internal and form- 
drag components challenge the best experimental tech- 
nician. The .results of such an investigation are 
greatly needed by air-frame designers to assist in the 
selection of the optimum air inlet from the overall 
standpoint. For example, although the recovery of 
a flush iniet is inferior to that of a full-protruding inlet, 
the reduction in form drag may, from the overall stand- 
point, offset recommended, 


therefore, that, wherever possible, all future air-inlet 


this inferiority. It is 


investigations include provisions for the determination 
of internal and form drag of the inlet. 

A word of caution should be given with respect to the 
use of the values of inlet-total-pressure recovery shown 
in Fig. 3 for internal-flow analyses. The inlet-total- 
pressure recovery presented on the curves is the inte- 
grated mean total pressure across the inlet (described 
under ‘‘Reduction of Data’). The total-pressure dis- 
tribution across the inlet varies, as may be seen in the 
typical plots in Fig. 8. It is well known that the effi- 
ciency of diffusers and bends depends fo a great extent 
on the inlet-total-pressure profile. This dependence, 


unfortunately, is not established for most duct geome 
tries but must await the results of systematic research) 
As a general design policy, however, it is believed ad¥ 
visable to precede any diffusion or bend with a length 
of constant-area duct in order that any irregularities jqy 
the velocity profile will tend toward symmetry. 


CONCLUSIONS AND RECOMMENDATIONS 


The results of a systematic investigation of the in 
fluence of design parameters on the performance of q 
These 
results are summarized in a form convenient for design 
purposes. It is concluded that investigations of thig 
type will help fill a gap in design information badly 
needed by aircraft designers. In so doing, much costly 
development testing of secondary air inlets may be 
avoided or minimized, and primary air-inlet configura- 
tions may be selected with a greater assurance of their 


certain class of subsonic air inlet are presented. 


success. 

It is recommended that the following research pro- 
grams be initiated without delay: 

(1) Systematic investigation of the influence of 
design parameters on the form drag of air inlets. 

(2) Systematic investigation of influence of inlet- 
total-pressure (or velocity) profiles on the performance 
of duct components such as diffusers, bends, ete. 

(3) Investigations of other basic inlet configura 
tions in a systematic manner, and the presentation of 
the results of these investigations in the form of design 
selection charts similar to those presented in this 
paper. 

(4) Extension of above to the high subsonic, tran- 


sonic, and the supersonic speed ranges. 
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